title: Adjusting Uniform Probability Distribution

text: I'm looking for a way to adjust the probability distribution of a uniform random function I'm using in a program. I want to find some discrete probability distribution that includes a parameter for "adjusting" how uniform the distribution is relative to the "center" event. For example, if that parameter, lets say _s=0_ , then the probability of 5 events is exactly equal at 20% each. If I were to change that parameter to _s=1_ , the function would now adjust the probability depending on the distance from the center. So, the probability would change to:   * Event 1: 16%   * Event 2: 18%   * Event 3: 20%   * Event 4: 22%   * Event 5: 24% And adjusting _s_ higher would yield increasing probability above the "center" event and decreasing probability below the center event, depending on distance from the center - all while retaining a total probability of 100%. Statistics is not my strong point, so I'm just hoping someone with some experience might understand what I'm looking for :)


title: Decision tree : handle attribute with many nominal values

text: I would like to build a decision tree from a training data. I have an attribute with many nominal values. For example, the department name attribute has about 20-30 values. I would like to group these values to 4-5 groups or any appropriate amount. How can I do this (in Weka is preferable)? Should I do this in preprocessing phase or/and learning phase?


title: How to select an error/weight for two-dimensional binned data?

text: Beginning with noisy data vectors $\mathbf{x}$ and $\mathbf{y}$, I have binned the data to vectors $\mathbf{x}_b$ and $\mathbf{y}_b$ of length $N_b$ with fixed linear ($\mathbf{x}_b^i - \mathbf{x}_b^{i-1} = \delta$) or logarithmic spacing ($\mathbf{x}_b^i/ \mathbf{x}_b^{i-1} = \rho$). I would now like to calculate a fit using, eg scipy.optimize.curve_fit. Prior to binning, the least-squares algorithm would look at all the points and values and would thus get decent weighting. However, now that I've destroyed that information, how can I produce good weighting/error factors for my data points? I saw a very thorough answer for a one-dimensional problem on this site but I don't know how to derive the two dimensional versions of the formulas. The difference is that the spread of points on the x-axis and the y-axis both need to be combined into one error parameter.


title: How to fit a simple count time series INAR(1) model in R

text: I am trying to perform a simple time series analysis with count time series data. My data is a sequence of small integer values like 0,1,2 and 3. I learned from various sources that INAR model would be appropriate with such data. My question is whether anyone knows R codes for fitting a simple INAR(1) model (regressing time series data on a binary dummy variable). Appreciate any assistance.


title: Independence in Linear Regression

text: This is probably a relatively simple question, but are b0 and b1 independent in a simple linear regression model?


title: how to measure the model recovery performance

text: I have some simulated linear model $y = X\beta+\epsilon$ where $\beta$ is sparse. I am comparing different techniques to recover the structure of $\beta$. I have so 4 values True Positives, False Positives, True Negative, False Negative. How can I summarize the performance with only one value?


title: How do I measure a speedup?

text: I have two versions of the same program. I apply a number of performance tests to both and measure their response time, for instance https://gist.github.com/valtih1978/d2cc2fe96fbbe1987ada               //  random(10000) //> compressed 10000-char msg into 66271 bits     //| in 59 (tree) vs 20 (table) msec     //| in 68 (tree) vs 17 (table) msec     //| in 45 (tree) vs 23 (table) msec     //| in 44 (tree) vs 19 (table) msec     //| in 42 (tree) vs 16 (table) msec     //| compressed 10000-char msg into 66319 bits     //| in 44 (tree) vs 15 (table) msec     //| in 77 (tree) vs 15 (table) msec     //| in 40 (tree) vs 95 (table) msec     //| in 47 (tree) vs 17 (table) msec     //| in 42 (tree) vs 15 (table) msec     //| compressed 10000-char msg into 66264 bits     //| in 45 (tree) vs 15 (table) msec     //| in 86 (tree) vs 15 (table) msec     //| in 39 (tree) vs 20 (table) msec     //| in 43 (tree) vs 14 (table) msec     //| in 41 (tree) vs 14 (table) msec     //| compressed 10000-char msg into 66315 bits     //| in 109 (tree) vs 17 (table) msec     //| in 85 (tree) vs 15 (table) msec     //| in 41 (tree) vs 21 (table) msec     //| in 94 (tree) vs 29 (table) msec     //| in 45 (tree) vs 30 (table) msec     //| compressed 10000-char     //| Output exceeds cutoff limit.      I want to figure out if one is faster and how much. This means that I have two samples, each has some variation. Now, how do I ensure that one is n times faster than another? Should I build some confidence distribution which will tell me that this it is k times faster with that probability? What can variance of rations can give me? Should I use resampling? I cannot even figure out how to compute average speedup. Should I take average of corresponding speedups, $speedup_{avg} = (b_1/a_1 + b_2/a_2 + \ldots)/n$ or somehow compute $b_1$ vs. all $a_n$s and then aggregate them? If I was interested in aboslute differences, I would take $speedup_{avg} = (b_1 - a_1 + b_2 - a_2 + \ldots)/n$ because it is equal to running all $a$ tests, then all $b$ tests and taking the avarage difference $speedup_{avg} = ((b_1 + b_2 + \ldots) - (a_1 + a_2 + \ldots))/n$ or average of resampling all $b_i$ against all $a_j$, $$speedup_{avg} = {1/n\ (b_1 - a_1 + b_1 - a_2 + \ldots) + (b_2 - a_1 + b_2 - a_2 + \ldots)/n + \sum_{i=3}^n{\sum{(b_i - a_j)}/n} \over n}.$$ The fact that averages of agerages converge seems not the case for the ratios, howver. For instance, $$avg_1 = {(b_1 + b_2)/(a_1 + a_2) \over 2} $$ $$avg_2 = {b_1/a_1 + b_2/a_2 \over 2} = {(b_1 + b_2 + b_2 a_1/a_2 + b_1 a_2/a_1) \over 2 (a_1+a_2)} = avg_1 + {b_2 {a_1 \over a_2} + b_1 {a_2 \over a_1} \over 2 (a_1+a_2)} $$ I am lost. Does it mean that any ratio-based statistic is not ergodic? I mean, can I just add up execution time of the first sample and divide by total executoin time of second to get the avarage ratio, can I just make one long run of the first program and long run of the second program instead of breaking them into pieces or I need that for the variance? Should I use harmonic measure (workloads are the same for all the tests)? Or should I use square root for averaging? How do I decide? Some people say that Fieller's Theorem must be used but I do not understand what should be proven here. Is there a right measure for my purpose or there are many? Which is the simples and does the job?


title: Should I categorise my continuous variable for use in binary logistic regression

text: DV: work response dummy, (1=household increased work, 0=hh did not increase work) IV: family size (continous variable) I am doing a logistic regression on these variables. My concern is that I know when using categorical predictors I need to make sure there are no empty cells. However, I'm not sure exactly how the logistic regression works with a continous predictor. I know that it will look at the change in the odds of the DV resulting from a one unit change in the continous predictor. However, does it treat the continuous variable in a similar way. Should I worry about empty cells in this case. For example, there may be only one family with 17 members. I hope my question is reasonably clear. I appreciate any help you could give me. Thank you.


title: Relative increase between two pairs of samples

text: Consider the following samples from four distributions: $a = \\{20, 2, 200\\} \qquad\qquad c = \\{1, 10, 100\\}$ $b = \\{22, 2.2, 220\\} \qquad\quad\:\,d = \\{1.2, 12, 120\\}$ I would like to say that the relative increase from $a$ to $b$ is significantly less than the relative increase from $c$ to $d$. * * * If the samples are assumed to be paired, with the $n$th element of $a$ corresponding to the $n$th element of $b$, this can be done with a Mann- Whitney U test on the element-wise ratios of the distributions: $\frac{b}{a} = \\{1.1, 1.1, 1.1\\}$ vs $\frac{d}{c} = \\{1.2, 1.2, 1.2\\}$ giving a p-value of 0.02 (tending towards 0 with larger samples). * * * But what if the samples are not assumed to be paired, or have different sizes? A 'ratio bootstrapping' test could be used instead: randomly select numbers $a_i$, $b_j$, $c_k$, $d_l$ with replacement from each distribution and see how often $b_j / a_i < d_l / c_k$. For this example, 31 of the 81 possible selections have $b_j / a_i > d_l / c_k$, giving a p-value of $\frac{31}{81} \approx 0.38$. Ratio bootstrapping is inconclusive here because the variance within each distribution is larger than the variability (e.g. difference in means or medians) between the pairs of distributions. Is there a statistical test which accomplishes the same thing as ratio bootstrapping, but is insensitive to intra-distribution variance and only focuses on the variability between distributions?


title: Calculating differences between groups and across time (age)

text: Let's say I have data of one variable $V$ measured for two groups of individuals ($a$ and $b$) and I want to see whether there is a significant difference between them and _over time_ , i.e. $V_a$ stays the same over time but $V_b$ decreases. Just to stress that $a$ and $b$ are two groups of individuals with different ages, and not two individuals measured over a period of time, how do I do that? I started reading things about ARIMA models but it gets quite technical, and I'm a fan of simplicity :-) Thanks much! omer


title: Single Categorical DV (3 levels), and Single Continuous Repeated Measure IV: which test?

text: I'm a Ph.D. psych student and am having trouble finding information on which test to use for a continuous repeated measures IV and categorical DV. I would really appreciate some help with this. The outcome variable is depressive symptom trajectories (classes derived from Latent Class Growth Analysis of depressive symptoms measured four times). There are three levels (classes) of this outcome variable. The predictor is cultural values, continuous data collected twice. I'd like to test whether a change in cultural values predicts which depressive symptom trajectory the participants belong to. Using Multinomial Logistic Regression, the results looking at only the first time-point of the IV are pretty interesting, and I'd like to build on this model by analyzing the IV as a repeated measure. It seems Repeated Measures Logistic Regression is only for repeated DVs though. I would be forevermore grateful for a nudge in the right direction. Kindly, Penelope


title: VaR in case of ARMA-GARCH?

text: How do I calculate the VaR in case of using an ARMA-GARCH approach? I am not good at time series, so I am more or less confused with the different possible notations of an ARMA-GARCH process. I hope the notation I use is correct. Suppose I use the following AR(1)-GARCH(1,1): \begin{align*} r_t=\delta + \phi_1 r_{t-1} + z_t\\\ z_t = \sigma_t \epsilon_t\\\ \epsilon_t \text{ is } N(0,1)\\\ \sigma_t^2=\alpha_0+\alpha_1 z_{t-1}^2+\beta \sigma_{t-1}^2 \end{align*} I mean, I would have to calculate a quantile at timepoint t. Normally in case of a more simpler approach: $r_t=\mu + \sigma_t*\epsilon_t$ and the $\epsilon_t$ is N(0,1) this would give $VaR=\mu + \hat{\sigma_t}*z_\alpha$ where $z_\alpha$ is the quantile of the standard normal. How to do this in case of the ARMA-GARCH approach? I now have no $\mu$, but an AR process? I read some papers, they all give more or less the ARMA-GARCH formula and say that they generate the forecast, but they give no formula, how to do this.


title: Bounded response variable [-1;1] - Should I transform it?

text: I am planning to use two response variables. One is bounded between 0 and 1, and I guess I can use a binomial (or related) error structure. The second variable is bounded between -1 and 1. I am not sure what should I do for this one. Should I transform it? Please note that I plan to use these variables in mixed-effect modeling (see related post on mixed-effect modeling: Mixed-effect modeling with paired observations & bounded response variables). Any help very much appreciated.


title: What current methods for statistical disclosure limitation are best trade-offs between data privacy and data utility?

text: In many situations raw microdata is not released by institutions due to privacy preserving. Many techniques are used for protecting sensitive values in data. But many of them can destroy multivariate structures in data and data itself is then not suitable for research at all. On the other hand, if data is not protected adequately, privacy preserving may not be assured. Therefore, any trade-off between privacy preserving and utility of data is desired. What methods are most popular for this purpose and attempt not to destroy multivariate structure of data while preserving privacy as much as possible?


title: Pros of Jeffries Matusita distance

text: According to some paper I am reading, Jeffries and Matusita distance is commonly used. But I couldn't find much information on it except for the formula below JMD(x,y)=$\sqrt[2]{\sum(\sqrt[2]{x_i}-\sqrt[2]{y_i})^2}$ It is similar to Euclidean distance except for the square root E(x,y)=$\sqrt[2]{\sum(x_i-y_i)^2}$ JM distance is claimed to be more reliable than Euclidean distance in term of classification. Can anyone explain why this difference make JM distance better?


title: Question about the probability chain rule

text: I've understood from this: Is this a correct statement of the probability chain rule? that in the chain rule for probability, conditioning can be done on different variables. I was wondering what sort of implications this has for language modeling where we try and assign a probability to a sequence of m words P(w_1, w_2, ..., w_m) or any sequence for that matter. For example, say we want to assign a probability to P(the, cat, ran). As seen here: http://en.wikipedia.org/wiki/Language_model#N-gram_models one would compute P(the, cat, ran) by calculating P(the)*P(cat|the)*P(ran|the,cat). But given that conditioning can be done on different variables, would it be correct to say that P(the, cat, ran) = P(ran)*P(the|ran)*P(cat|ran,the)? If so, how is the sequence of the words actually taken into account? The probability of "the cat ran" is obviously more likely than "ran the cat".


title: Machine learning with weighted / complex survey data

text: I have worked a lot with various nationally representative data. These data sources have a complex survey design, so the analysis requires the specification of stratification and weight variables. Among the data sources that are within my area of study, machine learning tools have not been applied to them. One obvious reason is that machine learning methods (currently) do not take into account weight and stratification variables. The goal of the weighted / stratified analyses is to obtain adjusted population estimates, which is different than the goal / purpose of machine learning. What thoughts do people have about using the nationally representative data sources and ignoring the weight and stratification variables? In other words, what would be your thoughts if you reviewing a machine learning study that was used nationally representative data but ignored the weight and stratification variables, assuming that the researcher / author was up-front about this methodological decision and was not making claims of nationally representative results? Thanks in advance!


title: Three-way interaction in multilevel model

text: I am doing a multilevel Analysis in which I test whether the interaction between two Level 1-Predictors (IV1 and IV2) is moderated by a Level 2 predictor (IV3). The Level 2 predictor is a dichotomous variable. The model I estimate reads as follows: Level-1 Model DVij = β0j + β1j*(IV1ij) + β2j*(IV2ij) + β3j*(IV1xIV2ij) + rij Level-2 Model β0j = γ00 + γ01*(IV3j) + u0j β1j = γ10 + γ11*(IV3) β2j = γ20 + γ21*(IV3j) β3j = γ30 + γ31*(IV3j) IV1 IV2 have been centered around the grand mean. I have two questions:   1. I find that the three-way interaction (γ31) is significant. I plotted a graph including four slopes. I wondered if there was a way to test which slopes significantly differ from one another, similarly to the approach of Dawson & Richter (2006) for multiple regression. Reference: Dawson, J. F., & Richter, A. W. (2006). Probing three-way interactions in moderated multiple regression: Development and application of a slope difference test. Journal of Applied Psychology, 91, 917-926   2. When I build the model, the interaction between IV1 and IV2 first is not significant (γ30). However, when I add the Level 2 predictor to the equation, the interaction between IV1 and IV2 becomes significant as well as the three-way interaction (γ31). How should I Interpret this pattern of results? Would my hypothesis that there is an interaction between IV1 and IV2 be supported or not?


title: How to compute sample size to detect a small change?

text: Suppose I wish to calculate the sample size needed to detect a one percent change in a proportion that is normally at 35%. I don't mean 35% to 36%, I mean 35% to 35.35%. Here is the calculation in R blowing up:               > library(pwr)     > p2 = 0.35     > p1 = p2*1.01     > p1     [1] 0.3535     > h = 2 * asin(sqrt(p1)) - 2 * asin(sqrt(p2))     > alpha = 0.05     > beta = 0.05     > pwr.2p.test(h=h, sig.level=alpha, power = 1-beta)     Error in uniroot(function(n) eval(p.body) - power, c(2 + 1e-10, 1e+05)) :        f() values at end points not of opposite sign      Here is a 3% calculation succeeding and saying the sample needs to be of size 53,991:               > library(pwr)     > p2 = 0.35     > p1 = p2*1.03     > p1     [1] 0.3605     > h = 2 * asin(sqrt(p1)) - 2 * asin(sqrt(p2))     > alpha = 0.05     > beta = 0.05     > pwr.2p.test(h=h, sig.level=alpha, power = 1-beta)               Difference of proportion power calculation for binomial distribution (arcsine transformation)                         h = 0.02194005                   n = 53990.98           sig.level = 0.05               power = 0.95         alternative = two.sided          NOTE: same sample sizes      I accept all the assumptions above (same sample size, two-sided test, etc.). See here and here for more info on pwr.


title: Determining appropriate sample rate for count data

text: **Context:** I am collecting Twitter data through a source that allows me to take a statistical sample of tweets that match a filter. I want to determine the appropriate sample rate in order to achieve a certain statistical confidence in my estimate of the actual number of tweets that would match that filter. **Sampler details:** For each tweet, the sampler calculates a new random floating point number uniformly distributed in the range [0,100], then compares that number to a given threshold (e.g. 1.6 for 1.6%). If the number is less than the threshold, the tweet is sampled; otherwise, it is discarded. **An example:** Let's say @justinbieber gets 500,150 mentions per day. I filter first to get all tweets that mention @justinbieber, then I sample at a probability of 1%. That should give me something like 5,001, which I'd then divide by 1% to get an estimate of 500,100 mentions/day. There's two sources of error here: 1) **Rounding error.** I am necessarily reducing the precision of my count data, then extrapolating back to the original. 2) **Sampling error.** While the sampler uses uniformly distributed probability in the range, I may not get exactly the percentage of tweets I am requesting. **The Goal** I want to calculate the appropriate sample probability I should use in order to achieve an estimate of twitter mention count with a certain confidence. For example, I might want to say with 95% confidence that my estimate is within 2% of the actual count. Given both sources of error here, what formulas or techniques allow me to calculate the proper probability to use in my sampler given an estimated per-day count? **Notes** The best solution should work equally well with small counts (<40) and large counts. It's certainly acceptable to sample at 100% for twitter handles with small counts, but there is a fixed cost per-tweet, and thus minimizing the number of tweets required is important when possible.


title: Simple regression models for data with a breakpoint

text: I am currently working on a segmented regression model with two variables $(x_i,y_i)$ for $i = 1..N$. The model should take the form: $y_i = \beta_0 + \beta_1 x_i$ for $x_i < x_{crit}$ $y_i = \alpha_0 + \beta_0 + \beta_1 x_i$ for $x_i \geq x_{crit}$ where $x_{crit}$ is a breakpoint value that has to be determined from the data, as well as the coefficients $\alpha_0, \beta_0$ and $\beta_1$. I have seen similar ideas in statistics / econometrics papers (i.e. segmented regression models / regression discontinuity design) but I would like to get some feedback on what kinds of model I should use. Ideally, I would like to use a well-documented framework where I can get some sort of confidence interval on the breakpoint $x_{crit}$.


title: In matlab, results from anova1 and multcompare disagree?

text: In matlab, I've ran               [p,stats,~] = anova1(X)      with some samples, and it indicates a significant difference in (at least) 2 means from the samples. i.e. $p < \alpha$ for some value of $\alpha$. However, running               multcompare(stats,'alpha',alpha)      doesn't indicate any significant difference between the data sets (i.e. there aren't any confidence intervals that don't contain 0). Is this a contradiction, or is this possible, and why? Many thanks.


title: Anomaly detection for one feature vector

text: I have a $n$-dimensional vector of ordered multiple testing $p$-values and I would like to reject the first values that are under a certain threshold $\alpha$. I am looking at this problem as an anomaly detection and clustering when there is an abrupt change that happens between rejected and non-rejected values. And so, I am wondering how to proceed to this problem when I have only one feature space which is my ordered p values. Do I have to map my vector onto a different feature space where I have bigger number of features or this is a simple bi-clustering problem using just one feature.


title: Logistic regression model fit assessment

text: I have run logistic regression on data. The concordant percentage is coming out 45. That is too low. However, Hosmer–Lemeshow test is coming out insignificant. It means we cannot reject the null hypothesis stating model fits the data well. That's good! The ctable also shows 70% of the data correctly predicted. How contradictory results of concordant and Hosmer-Lemeshow test should be considered? How should I interpret this model?


title: K-means Stuck in 1 cluster

text: I'm working on a problem using the encog Kmeans library and NO MATTER what features I add to the model, it always gets stuck in one of the clusters. ALL of the samples are lumped into one cluster every time. Does someone know if there's an error in the K-means code for Encog, or am I just doing something incredibly stupid? For the record, I've tried it with lots of different features that have different number ranges.


title: Possible inferences from a graph pattern

text: So, I had a weighted dynamic graph having info about 10 consecutive timesteps ( basically 10 files ). Now, I had to mine out patterns in the weight and structure of the complete graph. I did that. The output file (having around 10,000 rows) was something like,               Node 1 Node 2 Pattern      191    570    "00"        21     570    "00"       378    570    "00"       459    570    "00"       552    570    "00"       223    570    "00"      197    570    "00"       570    689    "00"       ...................      Basically gives 2 nodes, and the pattern associated. What I wish to ask is, what kind of model like (linear regression, or bar plots, .... ) can I use here, in order to gather some meaningful info, i.e. let's say, I am able to come up with a prediction model or I can say, that these are 2 graphs which are similar in nature.


title: How to test multiple group mean different than unconditonal mean using regression?

text: Let's say we have N apples from multiple vendors (A, B, C, D) and want to test if mean weight of each vendor group different than unconditional mean weight. I know I can use cross sectional dummy variable to test whether mean of B, C, D are different than mean of A, by creating 3 indicators Is B, Is C and Is D, making A base case. But here I want to use the population mean as the base case. Is there a way? I may need to add more factors to this model such as location and color Thanks


title: Choosing the number of clusters when only the distance matrix is available

text: I am performing a hierarchical clustering and a DBSCAN clustering on a distance matrix (A self-organizing map is something that I want to do as well). I am trying to find a method to choose the number of clusters based on some criteria (like an elbow plot with some criteria for its slope or a BIG criteria or…), but I am not sure what are the efficient and reliable methods of doing so when only the distance matrix exist. Any suggestion is appreciated.


title: error increasing with no of estimators in adaboost

text: My error gets increased when i increase the n_estimators value in               ada=AdaBoostClassifier(n_estimators=1, base_estimator=SVC(probability=True,kernel="rbf",gamma=.1,C=1))      for n_estmators the predictor is giving the least error, whereas as i increase the value of n_estmator then the error also increases what could be interpreted from this? thanks


title: Why multiplication and division are not allowed when using the interval scale?

text: **This is what I already know:**   There is no _true zero value_ on interval scales like temperature.   The zero value doesn't state that temperature becomes unavailable at that point.   Zero is just a human defined number which represents a certain level of temperature.   This is all because temperature is a relative scale. **What I'm confused about:**   * Why can't I say 40$^\circ\\!$C is twice as heat as 20$^\circ\\!$C?    * What about weight? It's relative, but I can say 40 kg weighs twice as much as 20 kg.


title: How to obtain the gamma distribution through convolution of two different distributions?

text: How can one obtain the gamma distribution through convolution of two different distributions? Could the gamma distribution be created as a non-trivial sum of $N$ random variables $X$ which have the same distribution and parameters? Trivial case is summation of fixed number $N$ of variables with gamma distribution as described on Wikipedia.


title: Variations of Assouad’s lemma

text: I'm reading the paper Optimal Rates of Convergence for Sparse Covariance Matrix Estimation by T. Cai and H. Zhou. They refer to a specific version of Assouad's Lemma (see lemma 2, p. 6 of the linked paper), which differs from the one I've found in Yu [1] and Le Cam [2]. I've also ordered van der Vaart [3], but it has yet to arrive. Can do you know of a resource that contains a proof of the version Cai states? Or can you give a proof? I'm especially interested because I'm having trouble with their proof of lemma 3 (p. 7 of linked paper, proof p. 21) and the lemmas seem very related. Specifically, I don't understand how the use of $\hat{\theta} = \mbox{argmin}_{\theta \in \Theta} d(T,\psi (\theta))$ in eq. (37) is meaningful as $T$ is a random variable. [1] B. Yu (1997). Assouad, Fano, and Le Cam. Festschrift for Lucien Le Cam. D. Pollard, E. Torgersen, and G. Yang (eds), pp. 423-435, Springer-Verlag. [2] Le Cam, L. (1986). Asymptotic Methods in Statistical Decision Theory. SpringerVerlag, New York. [3] van der Vaart, A.W. (1998). Asymptotic Statistics. Cambridge University Press, Cambridge.


title: Collapsing data for repeated measures ANOVA?

text: **Data** An experiment was run with 40 subjects. 20 subjects are in condition 1 and 20 in condition 2. Over the course of the experiment, subjects went through 60 trials. 20 trials with stimuli of type A, 20 with stimuli of type B and 20 trials with stimuli of type C. In each trial, 6 scores are measured (secondwise). **Image** ![](http://i.stack.imgur.com/WbZsY.jpg) Fig. Two plots, one for condition 1 and one for condition 2. Secondwise scores are plotted collapsed over subjects (20 per condition) but split according to trial type (A/B/C). Errorbars show SEM. **Are scores statistically different between condition 1 and condition 2?** Approach so far: Collapsed scores over type of trial per subject, so that there are 3*6 measurement points per subject instead of 3*20*6. Turned data into wide format and ran repeated measures ANOVA in SPSS. Only main effect of second is significant. Looking at the graphs, however, the blue one seems very different between the two conditions. **Question** Did I throw away information by collapsing data before running the ANOVA?


title: Dummy variable interaction

text: Consider the following model: $wage_i = \beta_0 + \beta_1male_i + \beta_2private_i +u_i$ where $male_i=1$ if male; $0$ if female, and $private_i=1$ if private school; $0$ if state school. The **base group is females from state schools**. So when $male_i=1$ and $private_i=1$, the coefficients $\beta_1 + \beta_2$ measure the average wage differential of males from private schools compared to females from state schools. But now suppose we add an **interaction between the two dummy variables** into the model so that: $wage_i = \beta_0 + \beta_1male_i + \beta_2private_i + \beta_3male_iprivate_i+u_i$ So now when $male_i=1$ and $private_i=1$, **what do the coefficients $\beta_1 + \beta_2 +\beta_3$ measure?** **What is the interaction between these two dummy variables adding to the model?** Thank-you!


title: MLE for mean of symmetric but otherwise unknown distribution

text: Given i.i.d. draws $x_1,...,x_n$ from $X$, where:   * $X$ has a finite mean $E[X]=\mu < \infty$,   * $X$ is symmetric about its mean, meaning $f_X(\mu+c)=f_X(\mu-c)$ for all $c$,   * The probability density function $f_X$ is not otherwise known.  Is it possible to prove the following? **Proposition.** The MLE for the mean of $X$ is the sample mean, $\hat \mu_{MLE}=\bar x = \sum_{i=1}^n x_i$. A proof or a counterexample would be great. I am willing to additionally assume that $X$ has a finite variance $Var[X]=\sigma^2 < \infty$, or other common basic assumptions, if that becomes necessary for the proposition to hold, or if it greatly simplifies the proof. I suspect that it may be possible to use the invariance of the MLE to transformations of the data to prove this, but it might follow from simpler facts about the sample mean.


title: What algorithm/model should I use for this kind of prediction?

text: I have many features with different kind of values, these features will lead to $(x,y)$. Example for features and their values:               F1 = Russia          F2 = 120 km          F3 = "Political instability"          F4 = 80%          F5 = "Winter"          F6 = 80000 ton      Features from F1 to F6 will produce (20, 75432). What model/algorithm is best fit to train the machine on features from F1 to F6 to predict $(x,y)$? And what is the best way to transform these different values to be crunched by the algorithm?


title: Example on Hidden Markov Model

text: I was studying Hidden Markov Model(HMM) recently. I was looking to cross check what I understood. I found a code on Forward-Backward and Viterbi is given in simple Python terms in Wikipedia. I could find some R("HMM"), Python ("NLTK", "Scikit-Learn"), Java ("JHMM", "Lingpipe") implementations. I got three questions here. (a) Wikipedia examples are nice but reviewing them I got a mixed response. Some sources are saying they are nice examples while some are doubting the implementations. How are they? (b) In NLTK HMM I tried to read source from the link below, I found it is given in http://www.nltk.org/_modules/nltk/tag/hmm.html. I found forward probability in def _forward_probability(self, unlabeled_sequence) and Viterbi algorithm in def _best_path(self, unlabeled_sequence) Could I identify them fine? (c) I am looking for a simple example of Forward Probability and Viterbi Algorithm implementation [preferably an example from text, great if I get a Python/R implementation also], --I could search the web but did not find much. May any one help me? If anyone of the esteemed members of the group may kindly help me. Apology for any wrong question. Regards, Subhabrata Banerjee.


title: Intepreting p-value of a regression coefficient

text: I was wondering how you interpret the p-value of a coefficient for a multivariate linear regression. Is the p-value intepreted only when everything else is fixed? Is there a way to incorporate the dependencies of other variables?


title: How to compare percentages for two categories from one sample?

text: I want to compare two percentages on a single categorical variable from one sample. For example, in my data set, a variable can take on two values (i.e. A and B). Of 586 instances, 87.4% is in Category A (512 of 586) and 12.6% is in Category B (74 of 586).   * Can I say 'my participants tend to be more concerned about Category A than B'?    * Is a statistical test required to make such a claim? If so What statistical test is required?


title: what if precision equals to 1?

text: I am trying to calculate precision for multi class classification for different datasets. but mostly I am getting precision value as 1. what does it mean? any wrong with my algorithm?


title: Simulate LMM data with specific error values using R

text: I have a need to simulate data in R using the following: Simulate data for 5 (or more) subjects, indexed by $i$, $(i=1,~...,~5)$ with $6$ observations per subject indexed by $t$ $(t=1,~2,~...,~6)$. Observations of the $Y$ variable are taken at various times $(\text{TIME}=2,~3,~5,~9,~11,~13)$, per subject using various assumptions about $u_{0i}$ , $u_{1i}$: $$ Y_{ti} = (β_0 + u_{0i} ) + (β_1+ u_{1i})*\text{TIME}_t + ε_{ti} \\\ \text{where }ε_{ti} \sim\mathcal N(0, σ^2),~σ^2=0.01,~β_0=5,~β_1=0.7 $$ This is a problem I have been given (it is for a non credit course that I am taking for fun). Don't want the "answer", but maybe an example that I can build on or a suggestion or two.


title: Do we need to report the median or the Mean when using a Kruskal-Wallis test?

text: I have two groups. There is a significant difference according to the KW test. But do I need to report the mean or median of the groups?


title: Weight variables to develop scoring algorithm

text: I am studying how to develop a prognostic index in medicine. I have read some articles and in all of them metholodgy is quite the same. Most of them use logistic regression, choose variables and build a model. Each variable in the model has a weight. For example, if the score has three variables (ie.fever, cough, pain) each of them adds a number of points and the total sum is the value of the prognostic index. If prognostic index is used to decide if a patient must be admited to a hospital you may say that fever present sums 2 points, cough 1 point and feel perfect substracts 1 point (adds -1 point). If prognostic index I superior than 2 patient is admitted and otherwise not. What I don´t understand after reading some articles using this methodology is how they weight variables. Do they use logistic regression obtained Odds to weight variable contribution to the index? Is it correct to give more weight to larger odds? Are there any other methods to do this? What kind of test would be used to test fitness of this prognostic model? Best regards,. Marcos


title: What is the difference between a frequentist approach with meta-analysis and a Bayesian approach?

text: Say I am performing an analysis looking at a particular health measure. I am interested in the difference in that measure between patients and controls and wether or not the difference is different from 0. There have been studies in the past looking at my same research question and health measure but in different samples of patients. In my Bayesian analysis I would build up a prior distribution based on the previous studies incorporating the mean difference and standard error. Please forgive me if this a newbie question as I am newly learning Bayesian stats, but **in what ways would the results from my Bayesian analysis be different from the results I would obtain from using an inverse variance weighted meta-analysis to combine the mean difference estimates from the prior studies with my current data**?


title: ROC curve for the same entries with different values

text: Maybe a silly question, but need your advice and comments. I have data entries (compounds) with two values, Value1 and Value2 (activity for two proteins). Is it meaningful to build the ROC curve for these data entries with class1 (zero) as value1 and class2 (one) as value2? So, if the ROC AUC is close to 0.5 - means that we cannot distinguish between the Class1 and Class2. Are those 2 proposals correct? Thank's a lot!


title: Find probability density intervals

text: I have the vector               x<-c(1,2,3,4,5,5,5,6,6,6,6,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,8,8,8,8,9,9,9,10)      (my actual vector has a length of >10.000) and i would like to find the intervals where the 90% of the density lies. Is the `quantile(x, probs=c(0.05,0.95), type=5)` the most appropriate or is there any other way? Thank you very much E.C


title: Limits of correction for over dispersion with Poisson / binomial data

text: Reading many help pages, I see that a common way to correct for overdispersion is to fit models using quasibinomial or betabinomial distributions. However, I can't find help concerning how far correction is achieved nor how to be confident that a reliable model has been fitted. I am interested in mutation data from many genes that fall into 10 ordered classes. Each observation comprises a number of bases surveyed (`N`) and the number of substitutions observed (`S`). Overdispersion is great (~22!) because the genes vary considerably in how fast they are evolving. There is also the issue that `N` also varies a lot from <100 up to >30,000. I want to test the hypothesis that average substitution rate (`S/N`) increases with class number. I have tried the following options:   1. (awful!) Ignore the fact the data are binomial, set a minimum N>=1000, and fit a Gaussian model with S/N as the response. Fit is poor, though arguably not dreadful.   2. Use cbind(S, N-S) and fit a model with a quasibinomial error structure. Seems to give a similar answer that avoids the crazily low P-values I get when using family=binomial, but how do I know if the degree of correction is valid? This seems a relatively simple problem, but I want to make sure the model is fitted properly and is, if anything, conservative.


title: Test for a logistic distribution in R

text: I have a set of data and I'd like to know whether this data set has a logistic distribution. When I made a histogram of my data set it seems to have a logistic distribution, but to be sure I'd like to test for a logistic distribution in R. ![enter image description here](http://i.stack.imgur.com/38kLD.png)   So my question is: **Is there a way to test your data for a logistic distribution in R and how do you do this?** Additional information: The data set consists of 8544 items. The data are horizontal distances in km between 2 geographical points. I need the data to construct a dispersal kernel for elephants. I could construct this dispersal kernel more easily if the data had a normal or a logistic distribution. I already used a Kolmogorov Smirnov test to test for a normal distribution (p<0.05), but the result was significant, so the data are not normally distributed.


title: Graphing effect size for coefficient of determination

text: If I have a significant correlation coefficient of r=.80 between Variable A and B, I can work out the effect size (coefficient of determination) by squaring it, which is 64%. I want to graph this in the simplest way possibe (given my non-statistical target audience). Can I use a 100% stacked bar graph for this purpose. This will show Variable B as 100% and on top of it would be Variable A which would be 64%. I can then graphically say that 64% of variance in Variable B can be attributed to Variable A. (Conversely, I can also graphically say that I am unsure of the remaining 36% of Variable B) The appeal of this approach for me is that I can show the effect size between A and B on a number of variables (e.g. gender, age, education) in one graph. This will also make a colourful presentation (which is good for a non statistical target audience!). I have seen some textbooks showing two circles, each representing a variable (e.g. A and B). The part of the circles that overlaps illustrates the effect size. I thought doing a 100% stacked bar graph was a better way. **From the discussion below, it appears that the scatterplot is the way to go on this matter. However, how do I show 64% on a scatterplot?** I think the point is being missed in the discussion below. It is easy to illustrate the relationship through the scatterplot but how is the effect size illustrated i.e. the actual perecentage as above. I can't see this percentage figure in any of the diagrams below.


title: Is there any difference between normal distribution theory and least square theory?

text: Different approaches are available for determining the effect of an IV(Independent variable) on another variable. How the two assumptions cited in the question affect the estimates.


title: Likelihood ratio test overestimates significance. More appropriate test?

text: I'm trying to find words that have a significantly distribution in one small portion of a large corpus than they do in the entire corpus. I'm using the standard likelihood-ratio test statistic $$\chi^2 = 2 \sum_{i} O_i log(O_i/E_i)$$ where the two classes are the small portion of the corpus and the rest of the corpus. The problem that I run into is that the small portion of the corpus is so much smaller than the whole corpus that the expected number of occurrences of any word in that small portion is on the order of 0.001. Therefore, the fact that a word occurs at all in the small portion is taken to be very significant, even if it's a very common word like a, the, of, or another stopword. Is there a way to do a better estimate that doesn't underestimate expected values so badly? Or would a different test be more appropriate here? Thanks so much for your help.


title: Poisson regression in a survival setting on a simulated set

text: I have a large dataset with patients and I'm studying a rare outcome (~ 2%) and death is a competing risk (mean age ~69 years). I've used the R "cmprsk" package for my statistics and it seems that competing risks and the Cox regression are performing similarly although the competing risk analysis is more conservative giving hazard ratios closer to 1. I've been suggested to do a Poisson regression on the data but the results don't make any sense and I would be really grateful to get some input on the benefits of doing this kind of analysis on survival data. I've created this simulation for creating a dataset with similar risk factors:               library("cmprsk")     # The time for the study     accrual_time <- 10     followup_time <- 1          base_risk <- list("event" = .015, "cmprsk" = .1)          risk_factors <- list(list("frequency"=.1,                      "event" = base_risk$event*.5,                      "cmprsk" = base_risk$cmprsk*2),             list("frequency"=.05,                      "event" = base_risk$event*1,                      "cmprsk" = base_risk$cmprsk*1),             list("frequency"=.05,                      "event" = base_risk$event*-.5,                      "cmprsk" = base_risk$cmprsk*0))          # Number of subjects     n <- 5000          # Create base time, sequential inclusion     time_in_study <- rep(c(1:n)/n*accrual_time + followup_time, 1)          set.seed(100)          # Create empty sets     x <- matrix(0, ncol=length(risk_factors), nrow=n)     time_2_event <- rep(0, n)     time_2_comprsk <- rep(0, n)          # Create each studied observation and outcome     for(i in 1:n){         # Set base risk         event_risk <- base_risk$event          comp_risk <- base_risk$cmprsk              for(j in 1:length(risk_factors)){             x[i, j] <- rbinom(1, 1, risk_factors[[j]]$frequency)[1]                  # If there is a risk factor defined             if (x[i, j] > 0){                 event_risk <- event_risk +                         risk_factors[[j]]$event                 comp_risk <- comp_risk +                          risk_factors[[j]]$cmprsk             }         }              # Time 2 event/risk is 1/rate meaning that higher number -> shorter time         time_2_event[i] <- rexp(1, rate=event_risk)[1]         time_2_comprsk[i] <- rexp(1, rate=comp_risk)[1]     }          cn <- c()     for(i in 1:length(risk_factors)){         ev_rsk <- risk_factors[[i]]$event/base_risk$event+1         cmp_rsk <- risk_factors[[i]]$cmprsk/base_risk$cmprsk+1         name <- paste("Risk factor no: ", i, "\n * ev=", ev_rsk, " cr=", cmp_rsk, " *", sep="")         cn <- c(cn, name)     }     colnames(x) <- cn          # Select the event that happens first: study ends, evenent occurs, a competing event occurs     time <- apply(cbind(time_in_study, time_2_event, time_2_comprsk), 1, min)          # Outcome identifiers     event <- (time_2_event == time) + 0     comprsk <- (time_2_comprsk == time) + 0     cens <- event+2*(event==0 & comprsk==1)          out.cox_ev <- coxph(Surv(time, event)~x)     summary(out.cox_ev)          out.crr_ev <- crr(time, cens, x, failcode=1)     summary(out.crr_ev)          out.cox_cmprsk <- coxph(Surv(time, comprsk)~x)     summary(out.cox_cmprsk)          out.crr_cmprsk <- crr(time, cens, x, failcode=2)     summary(out.crr_cmprsk)      The output makes sense but when I do a:               out.glm_pr <- glm(event ~ x, family="poisson")     summary(out.glm_pr)      It gives estimates of:   * RF 1 ~ .14    * RF 2 ~ .41    * RF 3 ~ -.23 My questions:   * Is the glm() code correct or should I somehow transform my data?   * Does the Poisson output make any sense and how should if so interpret it?   * What are the benefits/pitfalls in using Poisson regression for survival data? Thanks! * * * ## UPDATE After adding exp(out.glm_pr$coefficients) the results are almost identical to the competing risk regression, here's a forest plot that compares the three: ![A forestplot comparing the different methods - Poisson: 1.152  1.509  0.794, CRR: 1.151 1.524 0.812, Cox PH: 1.897 1.931 0.798](http://i.stack.imgur.com/14Zt0.png) The x-axis is perhaps not entirely valid (should be "incident rate ratios" for the Poisson regression) but why are the outcomes for CRR & poisson almost identical? As for testing over-dispersion I've found these two methods:               > library(qcc)     > qcc.overdispersion.test(event)          Overdispersion test Obs.Var/Theor.Var Statistic p-value            poisson data         0.9391878      4695 0.99902     >      > library(pscl)     > out.glm_nb <- glm.nb(event ~ x)     Warning messages:     1: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace >  :       iteration limit reached     2: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace >  :       iteration limit reached     > odTest(out.glm_nb)     Likelihood ratio test of H0: Poisson, as restricted NB model:     n.b., the distribution of the test-statistic under H0 is non-standard     e.g., see help(odTest) for details/references          Critical value of test statistic at the alpha= 0.05 level: 2.7055      Chi-Square Test Statistic =  -0.0139 p-value = 0.5       I conclude that there isn't any evidence of over-dispersion or are there other methods better suited for testing over-dispersion in this kind of survival data? The quasipoisson analysis gives similar values:               > out.glm_quasi_pr <- glm(event ~ x, family=quasipoisson(link="log"))     > round(exp(out.glm_quasi_pr$coefficients), 3)     (Intercept)       xRF 1       xRF 2       xRF 3            0.059       1.152       1.509       0.794


title: How to select from which machine came the request in a distributed service system?

text: Map/Reduce is great concept for sorting lots of data at once. What to do if you have small parts of data and you need to reduce it all the time? Simple example - choosing a service for request. Imagine we have `10` services. Each provides services host with sets of request `headers` and post/get arguments. Each service declares it has `30` unique keys - `10` per set.               service A:     name     id     ...      Now imagine we have a distributed services host. We have `200` machines with `10` services on each. Each service has `30` unique keys in their sets. But now to find to which service to map the incoming request we make our services post unique values that map to that sets. We can have up to or more than `10 000` such values sets on each machine per each service.               service A machine 1     name = Sam     id = 13245     ...     service A machine 1     name = Ben     id = 33232     ...     ...     service A machine 100     name = Ron     id = 777888     ...      So we get `200 * 10 * 30 * 30 * 10 000 == 18 000 000 000` and we get `500` requests per second on our gateway each containing `45` items `15` of which are just `noise`. And our task is to find a service for request (at least a machine it is running on). On all machines all over cluster for same services we have same rules. We can first select to which service came our request via rules filter `10 * 30`. and we will have `200 * 30 * 10 000 == 60 000 000`. So... `60` mil is definitely a problem... I hope to get on idea of mapping `30 * 10 000` onto some artificial neural network like Perceptron that outputs `1` if `30` words (some hashes from words) from the request are correct or if less than Perceptron should return 0. And I’ll send each such Perceptron for each service from each machine to gateway. So I would have a map `Perceptron <-> machine` for each service. Can any one tall me if my Perceptron idea is at least “sane”? Or normal people do it some other way? Or if there are better ANNs for such purposes?


title: How to compare effect of time across two groups on ten dependent variables?

text: I have one categorical independent variable with two levels (Group 1 and Group 2). These groups were measured at two times on 10 dependent variables, therefore, I have 2 sets of 10 variables for Time 1 and Time 2. For illustration sake, lets consider Hit_1 (tested at time 1) and Hit_2 (tested at time 2) These 10 variables are dependent, i.e. performance on one variable depends on performance on another variable (they are actually signal detection variables from a memory task ie hits, correct rejections etc). The 10 variables are all interval. Some of the variables were not normally distributed so I did a square-root transformation. How can I do the following:   1. test whether there is an effect for Group 1 (i.e. differences between Hit_1 and Hit_2)?   2. for Hit_1 test whether there is a significant difference between Group 1 and 2? For each of the above, how can I perform a multiple comparison analysis with an alpha correction?


title: Is Fisher's exact test right for this problem?

text: Let's say I have 4 plants and I measure `something` under 2 conditions. For example: The entries in the table show how many times I `saw something` under condition (C1 or C2) for plant (P1, P2, P3 or P4).                    P1,  P2,  P3,  P4      C1    0   20    0   19      C2  100   80  180  150       Now I would like to answer the question "what is the probability that these 39 times of seeing `something` is unevenly distributed between these 4 plants"? I am much interested in observation like this, rather than, say:                    P1,  P2,  P3,  P4      C1   25   20   15   19      C2  100   80  180  150       Is it right to use Fisher's exact test for this? I just extended the question from a 2*2 table to 2*4.


title: Asymptotic relative efficiency of median vs mean for Student t distribution

text: I am looking for a theoretical expression for the ARE of mean and median for student-t distribution against sample size (degree of freedom), exactly the blue curve shown in John Cook's blog. Can anyone point me to the right reference with the expression for it? Thanks


title: How does R's auto.arima() function determine the order of differencing when estimating a regression with seasonal ARIMA errors?

text: The authors' paper about the method (Hyndman and Khandakar 2008) states that an extended an extended Canova-Hansen test is used, although the helpfile updates this to say than an OCSB test is used. I have two questions, therefore:   1. What is an OCSB test?   2. Is the method of determining seasonal differencing the same for regression with ARIMA errors as it is for ARIMA modelling with no covariates?


title: Generalization and recall

text: I have some questions regarding generalization and recall:   1. It usually hear that a good classifier should "generalize" from the training data. Is that true?    2. I feel(?) that a "general" classifier should improve recall, however a "specific" classifier should have more precision, is this true?    3. If we think of a classifier that has to be learnt from a domain D1 and then used over a different domain D2, then a "generalized" classifier is preferred. In which other situations a "generalized" classifier is preferred?  Thanks in advance and sorry for the vague questions...


title: The ethics of using an optimal multiclass feature set for binary classification

text: I'm currently trying to find the best feature set/network architecture configuration for a binary classification problem, however to approach it via the usual means of building and testing does not seem feasible. Since I need to test four possible network shapes (with differing layers and sizes) over two different types of network (BLSTM and MLP), where there are six different feature sets each with 36 different classes, it quickly becomes quite messy, resulting in 2*4*6*36 = 1728 possible networks. I've decided to cut out the factor of 36 by training a bunch of networks capable of multiple classification, and then choosing the one which returns the best results to be the template for my binary classifier. I'm wondering though, is this a valid approach? Is it reasonable for me to assume that the network/feature set which gives me optimal performance on a multi-class problem will also work for the binary case? If so (or no), would anyone be able to point me in the direction of any helpful papers? Thank you!


title: Partial correlations among three distance matrices

text: I am trying to test the hypothesis that related species should deviate more in niche space the more they overlap in geographic space. In other words, related species either forage similarly and don't show a lot of geographic overlap, or if they overlap a lot geographically, then they have shifted a lot in niche space. Phylogenetic distance (how closely the species are related) explains a good deal of the variance in a distance matrix of how far apart the species are in niche space. So, I have three distance matrices among species: genetic distance, distance in niche space, and average geographic overlap. I don't think a partial Mantel test is exactly what I want, but I could definitely be wrong. Semi- mathematically, I'm thinking there is a strong correlation between the genetic distance matrix and distance in niche space (i.e. closely related species are not far apart in niche space), BUT when this doesn't hold true, it's because there is a large value in the corresponding cell in the geographic overlap distance matrix. What is the correct statistical approach to employ? Any specifics about how to run this in R are welcome but not required.


title: Choosing between two parameters in a model

text: I have a few parameters that are related (let's call them X1 and X2), and I want to use whichever one will provide the strongest model. The model has many other parameters. Would I simply be able to compare the AICc of these two models? Model using X1:               Y ~ X1 + X3 + X4 + X5 + X6       Model using X2:               Y ~ X2 + X3 + X4 + X5 + X6      This is confusing me as I'm not sure what other parameters will ultimately belong in the "correct" model. Ultimately, X3, X4, etc. may be thrown out. So would I want to test this with potentially "bad" parameters?


title: Perform analysis about model prediction diferrences using CARET

text: Using **CARET** package one can perform an analysis on differences between various models obtained using a dataset (the training dataset, trainSet) and the model that best fits the trainSet: For example, two models were determined using _train_ function, marsFit and rfFit:               resamps <- resamples(list(MARSCV = marsFit,RFCV = rfFit))     summary(resamps)     #Since models are fit on the same versions of the training data we can compute the differences, then use a simple t-test to evaluate the null hypothesis that there is no difference between models.     modelDifferences <- diff(resamps)     summary(modelDifferences)     bwplot(modelDifferences, layout = c(2, 1),            scales = list(x = list(relation="free")))      Is there a way to do a similar analysis (t-test to compare mean values for example) using the model predictions on the test set (separate dataset than trainSet)? **CARET** has a _predict_ function from which one can evaluate the differences between the observed (testSet) and the predicted (model output) values. Thanks!


title: A simple application of conditional probability

text: I am practicing using conditional probability and have the following question > Suppose the probability of contracting a rare tropical disease Chickungunya > on a holiday to India is 0.0001. Fred, on his return to the United Kingdom > from a Jungle Safari in India, feels unwell and decides to see his doctor. > Fred complains to his doctor of joint pain. The doctor had heard that there > had been some cases of Chickungunya reported in India and also knew that 64% > of people with Chickungunya complained of joint pain. There are other causes > of Joint Pain however and the doctor estimates that the probability of joint > pain in the absence of Chickungunya is 0.6. The doctor ordered a very > reliable test used to diagnose Chickungunya that returned a positive result > 99% of the time when the patient had Chickungunya and only returned a > positive result 4% of the time when the patient tested was not suffering > from the disease. The test result returned for Fred was positive. Then I am asked three questions: > (a) Based on the evidence to date can the doctor conclude that Fred has > Chickungunya? > > (b) Realising that there are no other symptoms of the disease or any other > more accurate tests available, the doctor decides to ask for the test to be > repeated. Once again a positive test result is obtained. How, if at all, > will this affect his belief in the diagnosis of Chickungunya? > > (c) Can the strategy of repeating the same test over and over again improve > his belief in a positive diagnosis of Chickengunya? If so, what is the > minimum number of times that the test must be carried out, assuming a > positive result each time, for the belief in Chickungunya to be greater than > the belief in not Chickungunya? Here is my attempted solution for part (a): First I extracted all the information out of the given text using the notation $P(d),P(JP),P(+)$ as probability of disease, probability of joint pain and probability of getting a positive on the test respectively. $$P(d)=0.0001\quad P(JP|d) = 0.64$$ $$P(JP|¬d) = 0.6 \quad P(+|d) = 0.99$$ $$P(+|¬d) = 0.04$$ $$P(d|(+, JP)) = \frac{P((+, JP)|d)P(d)}{P(+, JP)} = \frac{P(+|d)P(JP|d)P(d)}{P(+,JP)}$$ I have all the numbers above apart from $P(+, JP)$ but I figured you could calculate it by looking at disjoint cases of having the disease or not and get $$P(+, JP) = P(d)P(+|d)P(JP|d) + P(¬d)P(+|¬d)P(JP|¬d)$$ Is this allowed? Then the answer comes out as 0.0026. For part (b) I know intuitively what is happening but I can't put it into symbols can anyone give me a suggestion? I think if I can solve part (b) I should be able to work out (c).


title: Error: object 'descr' not found

text: I am trying to follow up the example code in the "Building Predictive Models in R Using the caret Package" paper from Max Kuhn[1]. Here is the part of the code:               library("caret")     set.seed(1)     inTrain <- createDataPartition(mutagen, p = 3/4, list = FALSE)      However, after the following,               trainDescr<-descr[ inTrain, ]          Error: object 'descr' not found      I get the "Error: object 'descr' not found" message. Looks like 'descr' is not recognized. Can anyone tell me how to fix this please? [1]: Kuhn, M. (2008),   "Building Predictive Models in R Using the caret Package,"   _Journal of Statistical Software_ , November, Volume **28** , Issue 5.   link


title: Equivalence in model performance metric like RMSE

text: I calculated the root mean squared error (RMSE) to compare the simulated values of a hydrologic model with the corresponding observations for three observational datasets.The number of data points in these datasets are 30, 40, and 60. How to find out if the RMSE values calculated with respect to these datasets are equivalent ? Does the change in the size of the datasets affect the power of the performance metric? How do I detect if it does? Is there a testing mechanism? Any suggestions ? Thanks.


title: Price levels and inflation in regressions

text: How would we include price levels in a regression? Would we just take the log of the price index number, e.g.: Consumer Price Index number provided by the Office of National Statistics (UK statistical agency)? If not, then how would we include it? Is it acceptable to just put in the log of the inflation rate as a regressor on the RHS? If I am modelling savings and want to see how it responds to price levels, how would I do this?


title: Two groups measuring their blood conc... significance between two group and time?

text: hello I am doing analysis on a research and need to do the following We have two groups (control and patient) we are measuring their blood concentration every 10 minutes for 48 times (repeated measurement). What stats analysis should we use. Seeing if the the two group( either control /patient) have a significance in the blood concentration and also if there is a significance/relation between the blood concentration measurements (fist 10 vs 20 vs ~480min) Iam guessing we need to use Generalized linear equation or MANOVA ... or other?. Thank You


title: How to sum two variables that are on different scales?

text: If I have two variables following two different distributions and having different standard deviations... How do I need to transform two variables so that when I sum the two result is not "driven" by more volatile one. For example... Variable A is less volatile than variable B (ranges from 0 to 3000) and variable B goes fro. 300 to 350. If simply add the two variables together the result will obviously be driven by A.


title: Simple linear regression output interpretation

text: I have run a simple linear regression of the natural log of 2 variables to determine if they correlate. My output is this:               R^2 = 0.0893          slope = 0.851          p < 0.001      I am confused. Looking at the $R^2$ value, I would say that the two variables are **not** correlated, since it is so close to $0$. However, the slope of the regression line is almost $1$ (despite looking as though it's almost horizontal in the plot), and the p-value indicates that the regression is highly significant. Does this mean that the two variables **are** highly correlated? If so, what does the $R^2$ value indicate? I should add that the Durbin-Watson statistic was tested in my software, and did not reject the null hypothesis (it equalled $1.357$). I thought that this tested for independence between the $2$ variables. In this case, I would expect the variables to be dependent, since they are $2$ measurements of an individual bird. I'm doing this regression as part of a published method to determine body condition of an individual, so I assumed that using a regression in this way made sense. However, given these outputs, I'm thinking that maybe for these birds, this method isn't suitable. Does this seem a reasonable conclusion?


title: Unequal sample sizes for two-way repeated measures

text: I have 3 different treatments with the observations of 6 different repetitions. However, all these different treatments (groups) have different numbers of patients. The first group has 26 patients, the second group has 28 patients and the third group has 19 patients. Thus, I can say that I have unequal sample sizes for two-way repeated measures anova. When I do the analysis by using SPSS, it calculates the sum of squares and degrees of freedom by using the minimum sample size of the third group, which is 19. Solving the problem and doing the comparisons by hand is time consuming. I would like to know how to handle this kind of problem, which is unequal sample sizes in each groups, by using the computer. Is there any other program, (SAS or R?) that deals with that, apart from SPSS, which obviously can not?


title: Non-overlapping state and measurement covariances in Kalman Filter

text: I have implemented an EKF to track the state of a moving vehicle in 3D. My question is twofold:   1. What happens if a measurement (e.g., for $y$ velocity) has a covariance that does not overlap with the current state and covariance? As an example, the measurement may have a value of 4 with an variance of 0.2, but the state vector's value for $y$ is 6 with a variance of 0.3.    2. The reason I ask question 1 is that I am getting negative values along the diagonal in my covariance matrix, signaling negative variances, which is clearly wrong. I have tried using the Joseph form update equation, i.e.,  $ P' = (I - KH)P(I-KH)^T + KRK^T$ However, my covariance matrix is showing negative values after the predict stage, not the update. The issue is with the projection of the covariance matrix using $P' = APA^T + Q$ where $A$ is actually a Jacobian matrix of my transfer function. Interestingly, the Wikipedia article has this for an EKF: $P' = AP + PA^T + Q$ I have raw debug output matrices that I can post if needed. The issue seems to be some negative covariance values between my linear and rotational velocities. Apologies in advance if this is a simple question. I'm familiar with Kalman Filters, but have never implemented an EKF that was this complex before.


title: Posterior conditional probability derived from Bayesian learning

text: Please refer to my previous question for the context. Now a new patient with disease class $Y'$ has measured severities $x'$ is given. What is the posterior conditional probability that $Y'=1$ given $x'$ and the original data $D$? My approach is $$ \begin{align*} \mathbb{P}(Y=1\mid x',D)&=\int\mathbb{P}(Y=1\mid x',\pi,\theta,D)\mathbb{P}(\pi,\theta\mid D) d\pi d\theta\\\ &=\int\mathbb{P}(Y=1\mid x',\pi,\theta)\mathbb{P}(\pi\mid D)\mathbb{P}(\theta\mid D)d\pi d\theta \\\ &=\int\frac{\mathbb{P}((X=x',Y=1)\mid\pi,\theta)}{\mathbb{P}(X=x')}\mathbb{P}(\pi\mid D)\mathbb{P}(\theta\mid D)d\pi d\theta \end{align*} $$ But then it seems dead. Could anyone help? BTW, the mean of Beta distribution and Dirichlet distribution is given in the question, I thought they might be useful.


title: GRS statistic in R to test that intercepts from multivariate panel regression are jointly zero?

text: The GRS statistic is the Gibbons et al. (1989) statistic that tests whether the estimated intercepts from a multiple regression model are jointly zero. The typical scenario involves a multivariate linear panel regression where you are explaining the returns to securities in terms of its exposures to factor return series. Theoretically, a good factor model will have an intercept statistically indistinguishable from zero. How do I calculate the GRS statistic in R? Thank you


title: Issues using Hermite approximation to bivariate Gaussian in R

text: In _"On Gaussian-like Densities of Order Greater than Two"_ (Willett, P. Thomas, J. B., 1987), section II, the author state: $\mathcal{N}(x,y,\rho)=\phi(x)\phi(y)\sum_{n=0}^{\infty}\rho^nH_n(x)H_n(y)$ where   1. $\phi(.)$ is the unit normal density,   2. $H_n(.)$ is the $n^{th}$ Hermite polynomial,   3. $\mathcal{N}(x,y,\rho)$ is the bivariate Gaussian density with correlation $\rho$. I tried to repeat this (approximately, that is up to the 10th Hermite polynomial) in R:               library(PolynomF)     x<-polynom()     H<-polylist(1,x);for(n in 2:10) H[[n+1]]<-x*H[[n]]-(n-1)*H[[n-1]]     Hp<-as.function(H)     #The Hermite polynomial of order 10     rho<-0.7     R<-c();for(n in 0:10) R[[n+1]]<-rho^n          HPA<-function(z) prod(dnorm(z))*sum(Hp(z[1])*Hp(z[2])*R)          z<-runif(2)     HPA(z)      Which is completely off the mark (i.e. for $z=(0.65,0.63)$ i get 1194). What am i doing wrong ?


title: How do I use the SVD in collaborative filtering?

text: I'm a bit confused with how the SVD is used in collaborative filtering. Suppose I have a social graph, and I build an adjacency matrix from the edges, then take an SVD (let's forget about regularization, learning rates, sparsity optimizations, etc), how do I use this SVD to improve my recommendations? Suppose my social graph corresponded to instagram, and I was tasked with the responsibility of recommending users in the service, based only on the social graph. I would first build an adjacency matrix A (mxm), take the SVD, A = U s V, choose the first k eigenvalues, then what? I would presumably create a new set of matrices: U_new ~ mxk s_new ~ kxk V_new ~ kxm then what does one do? I've looked on the web, and most links focus on calculating the SVD, but no one tells you what to do with it. So what should I do?


title: Categorical Dependent Variable in R

text: How does one interpret the beta coefficients in R for a model that has a categorical dependent variable that is not binary? For instance, in the hypothetical model below, I would interpret B1 as being both statistically significant and positive, but if my categorical dependent variable can take on 4 categories, then which category does B1 positively predict? That is, if the DV can be (blue, red, green, yellow), what does B1 predict? As an side, I am dealing with a very small sample size.                  bayesglm(formula = DV ~ B1 + B2 + B3, family = binomial(link = "probit"),          data = data)          Deviance Residuals:           Min        1Q    Median        3Q       Max       -1.07714   0.03221   0.14180   0.25680   0.86569            Coefficients:                                  Estimate Std. Error z value Pr(>|z|)       (Intercept)                -6.336e-01  1.859e+00  -0.341   0.7332       B1                          2.691e+00  1.347e+00   1.998   0.0457 *     B2                          1.586e+00  1.297e+00   1.223   0.2212       B3                          6.637e-01  1.755e+00   0.378   0.7052


title: Comparing facial curvature of two human faces

text: I need help determining if what I'm doing is sufficient. I have extracted the "curvature" from two human faces (essentially the curve that goes from one ear, down the face to the chin, and back up to the other ear). Now I want to compute a curve similarity with another human face. To try to make this invariant to the possible different curve sizes (since the images these faces are taken from can be taken from varying distances from the person) I've normalized both the X and Y dimensions by subtracting the mean and dividing by the standard deviation (I'm not sure this is correct). Finally, I'm using Dynamic Time Warping to compare these curves, this is the kind of result I get: ![DTW on Face Curves](http://i.stack.imgur.com/Ef82s.jpg) Basically I'm asking whether this approach sounds reasonable.


title: How to generate uniformly distributed points in the 3-d unit ball?

text: I have posted a previous question, this is related but I think it is better to start another thread. This time, I am wondering how to generate uniformly distributed points inside the 3-d unit sphere and how to check the distribution visually and statistically too? I don't see the strategies posted there directly transferable to this situation. Please help. Many thanks.


title: Calculating the marginals given the conditional distributions

text: Suppose $X$ and $Y$ have conditional distributions given by: \begin{align} f(x|y)&\propto ye^{-yx}\;\;\text{for}\;\;0<x<B<\infty\\\ g(y|x)&\propto xe^{-xy}\;\;\text{for}\;\;0<y<B<\infty \end{align} According to the article, the marginal distribution $g(x)$ is not easy to calculate, but it exists since $B<\infty$. It is easy to determine that: \begin{align} f(x|y)&= \frac{e^B}{e^B-1}ye^{-yx}\;\;\text{for}\;\;0<x<B<\infty\\\ g(y|x)&= \frac{e^B}{e^B-1}xe^{-xy}\;\;\text{for}\;\;0<y<B<\infty \end{align} Now since $f(x|y)f(y)=f(x,y)=g(y|x)g(x)$, this means that: $$\frac{f(y)}{g(x)}=\frac{g(y|x)}{f(x|y)}=\frac{x}{y}.$$ And thus that: $$\frac{1}{g(x)}=\frac{1}{g(x)}\int_0^Bf(y)dy=\int_0^B\frac{f(y)}{g(x)}dy=\int_0^B\frac{x}{y}dy=x\ln(y)\Big|_{y=0}^{y=B}=-\infty.$$    This would imply that $g(x)=0$, and by a similar argument, that $f(y)=0$. I'm not entirely sure what to make of this. Am I to conclude that since the ratio of marginals is an indeterminate form (zero over zero), that this method of calculating them fails in this particular case? and that they aren't in fact both zero. Also, is this a common technique for calculating marginals given the conditionals? I came up with it while trying to see why the technique of Gibbs sampling couldn't be replaced by the (possibly numerical) evaluation of some (possibly quite difficult) integral. Finally, like I mentioned above, the article says the marginals are not easy to calculate, which seems to imply that one could in theory calculate them, how would this be done?


title: Determining the amount of lag in an autoregressive model

text: I have done a lot of work in regression (time-invariant) but I am just now studying forecasting. My question is about determining the amount of lag to use in an autoregressive model. I assume that in forecasting you use your data up to time t to train your model and then do an out-of-sample test on the data after time t to test your model's performance. Is there such a thing as overfitting in autoregressive models? That is, using too much lag such that the model does not generalize well to out-of-sample data. Or is it that you stop seeing an improvement in performance after a certain number of lags?


title: R - Robust Ancova: Documentation/Examples for the ancmg function in the WRS package

text: Could anyone point me to a documentation or even better an example for the ancmg1 or ancmg functions of the WRS-package. There is a very short paragraph in the Wilcox' book about Robust Estimation, which describes very briefly the parameters that have to be supplied, but nothing about the actual method used or assumptions that have to be satisfied. In particular, I would like to know what the error message               simpleError in near(xm, pts[i], fr = fr): All measures of dispersion are equal to 0      means if I perform a robust Ancova with the function ancmg1 for a single covariate.


title: cross-validation and standardization

text: In this thread @whuber gave detailed answer about using training data statistics for standardizing cv dataset. My question is how to standardize hold-out dataset in n-fold cross-validation if some features are polynomial? Standardizing polynomial features will render them useless, right? Thanks


title: How do I analyze a case of breaking a cable at any point given that the breakage area length is zero?

text: Suppose I have a long network cable buried into ground and construction workers will sometimes unintentionally damage it when doing excavation works. I need to analyze the probability of breakage assuming it happens randomly. If I try to address the probability of breaking the cable at any specific point then it will be zero because the "length of a point" is zero so it's kind of impossible to unintentionally damage the cable at a specific point. So I have to divide the cable into sections for the analysis so that I then have sections of non-zero lengths and so non-zero probabilities of damage. How do I choose lengths of sections and what other ways to address the probability calculation are there that deal with "length of breakage point" being zero?


title: Learning from relational data

text: **Settings** Many algorithms operate on a single relation or table, while many real-world databases store information in multiple tables (Domingos, 2003). **Question** What types of algorithms learn well from multiple (relational) tables. In particular, I am interested in the algorithms that are applicable to the regression and classification tasks (not the network analysis oriented ones, e.g. link prediction). * * * I am aware of several approaches listed bellow (but am sure that I am missing some):   * Multi-Relational Data Mining (MRDM) (Dzeroski, 2002)   * Inductive Logic Programming (ILP) (Muggleton, 1992)   * Statistical Relational Learning (SRL) (Getoor, 2007) * * * Džeroski, S. (2003). Multi-relational data mining: an introduction. ACM SIGKDD Explorations Newsletter. Getoor, Lise, and Ben Taskar, eds. Introduction to statistical relational learning. MIT press, 2007. S. Muggleton and C. Feng. Efficient induction of logic programs. In Proceedings of the First Conference on Algorithmic Learning Theory, pages 368–381. Ohmsha, Tokyo, 1990.


title: Assumption tests in discrete-time hazard model

text: I am currently working with a random coefficient discrete time hazard model. The following assumptions are subject to this sort of model:   * linear additivity assumption   * proportionality assumption   * no observed heterogeneity assumption Now I would like to test these assumptions somehow. In the papers I have read so far, these assumptions were always taken as given, and I have not yet seen any standardized tests for them. Do any standardised tests exist for these assumptions or is rather a matter of the right specification of the model by trying out and comparing the deviance etc. of the different model specifications?


title: Degrees of freedom in a simultaneous confirmatory factor analysis

text: I am performing a simultaneous confirmatory factor analysis to check whether the items of a survey refer to intended separate underlying concepts without cross loading on the concepts. Since I am doing this in the context of a course on SEM, I am also paying attention to how this is estimated by the software. I learned that the model should not be under-identified, where there are more parameters to estimate than there is information available, i.e. the model is more complex than the data. Now what I think is, that my model is under-identified, however, the `lavaan` package in `R` has no problem fitting it, and suggests there are even 38 more degrees of freedom. My model looks like this: ![enter image description here](http://i.stack.imgur.com/3Jpmf.png) The variables with an asterisk are ordinal manifest variables which are turned into continuous latent response variables. This is the model specification in R:               TP.cfa <- 'tpeff =~ plcpvcr + plccbrg + plcarcr                  tpfairdis =~  plcvcrp + plcvcrc                 tpfairproc =~  plcrspc + plcfrdc + plcexdc                 cooperation =~ caplcst + widprsn + wevdct'     fit.cfa <- cfa(TP.cfa, data=ESS5BE, ordered=c("plcvcrp","plcvcrc", "plcrspc", "plcfrdc", "plcexdc", "caplcst", "widprsn", "wevdct"))     summary(fit.cfa, fit.measures=TRUE, modin=T)      If I count correctly there are 25 pieces of information available: 3 for the 3 continuous manifest variables {k(k+1)/2}, and 22 proportions for the ordinal manifest variables (6 of these have 4 response categories, 2 have 3 categories, and for one ordinal variable c-1 proportions are unique information). Lavaan estimates 53 parameters: 6 covariances, 7 factor loadings, 3 intercepts, 22 thresholds and 15 variances. Now my question is, would R be able to estimate the model, even if it is under-identified? Or do I really have an over-identified model, with 38 degrees of freedom, and then where do these come from?


title: demonstration of benefits of ridge regression over ordinary regression

text: I can understand **ridge regression** is better than **ordinary regression** in case of **multiple-collinarity**.               set.seed(123)      x1 <- rnorm(100)      set.seed(223)      x2 <- rnorm(100,mean=x1,sd=.01)      set.seed(344)      y <- rnorm(100,mean=3+x1+x2)      myd <- data.frame(x1, x2,y)      cor(myd)               x1        x2         y     x1 1.0000000 0.9999383 0.8701704     x2 0.9999383 1.0000000 0.8702772     y  0.8701704 0.8702772 1.0000000      The ordinary regression:               lm(y~x1+x2)$coef     (Intercept)          x1          x2         3.063971   -1.027534    3.097806       The ridge regression using `MASS` package               require(MASS)     lm.ridge(y~x1+x2,lambda=1)                    x1       x2      3.066382 1.019101 1.043206       The lambda (shrinkage) = 0 is equal to ordinary regression.                lm.ridge(y~x1+x2,lambda=0)                      x1        x2       3.063971 -1.027534  3.097806      **QUESTIONS:** (1) How can we show that **ridge regression is performing better than ordinary**? any estimators like error ? How to calculate ? (2) Are there other situations where ridge is better OLR ?


title: Exponential distribution qq plot homework question

text: For an exponential qq plot, we fix the theoretical distribution to have λ=1. In order to compare the plotted points with the line y=x, what kind of rescaling should be applied to the data and why? “Hint“ given with this problem: If X follows an exponential distribution with parameter λ, then λX follows an exponential distribution with parameter 1. I don’t know how to use this hint in answering this question. This is a homework problem. So far, I have considered the graph of ln (x) but not sure what to do with that and the hint.


title: Obtaining posteriors for exclusive and exhaustive hypotheses

text: I'm trying to solve a Bayesian problem where I have two mutually exclusive and exhaustive hypotheses: $H_1$ and $H_2$. Given Baye's formula: $$P(H|D) = \frac{P(H)P(D|H)}{P(D)}$$ (where $D$ is my data) I've read that one way to obtain the posteriors $P(H_1|D)$ and $P(H_2|D)$, is to create a table like this one:                     P(H) (prior) | P(D|H) (Likelihood) |         ____________________________________________________     H1 |      P(H1)    |        L(H1)        | P(H1)*L(H1)     H2 |      P(H2)    |        L(H2)        | P(H2)*L(H2)      Then summing the values in the third column I get the _normalization constant_ $P(D)$ as: $$P(D) = P(H_1)L(H_1) + P(H_2)L(H_2)$$ which yields for each hypothesis the posteriors as: $$P(H_1|D) = \frac{P(H_1)L(H_1)}{P(H_1)L(H_1) + P(H_2)L(H_2)}$$ $$P(H_2|D) = \frac{P(H_2)L(H_2)}{P(H_1)L(H_1) + P(H_2)L(H_2)}$$ My knowledge in statistics is not much and I haven't been able to confirm this anywhere else. Could someone tell me if this is correct?


title: What kinds of PRNG exist outside Linear congruential generators?

text: I'm kinda running out of words, Linear congruential generator are really low quality PRNG but it was all I needed for some basic stuff, now I need to fill the blank and get to know other families of PRNG but I don't have any idea about what are all the available families of PRNGs. Someone can possibly write the hierarchy or the list of all the possible kinds of PRNGs out there ?


title: What is the distribution of a chi-square minimizing function?

text: Suppose I have a set of $N$ experimental points of the form \begin{equation} \\{x_i, y_i, d_i\\}, \end{equation} where $i=1,...,N,$ and $d_i$ are errorbars for $y_i$. To fit the data, I minimize the reduced chi-square \begin{equation} \chi^2(p) = \sum_{i=1}^N \frac{[y_i - f(x_i,p)]^2}{d_i^2}, \end{equation} where $f(x,p)$ is a (generally non-linear) function parametrized by some parameter $p$ (there might be more than one parameter, but it doesn't really matter). My question is: given the optimal parameter $p_0$, i.e. $\chi^2(p)$ is minimal at $p=p_0$, and assuming the $y_i$'s are independent and are Normally distributed, what can be said about the distribution of $f(x, p0)$?


title: Using lead-lag relationships for time series prediction

text: I have a large number (hundreds to thousands) of noisy time series that represent contemporaneous observations from different subjects. I hypothesise that there exist lead-lag relationships between observations for different subjects (or groups of subjects.) I would like to explore the potential use of such lead-lag relationships for the purposes of predicting the individual series. What methods might I consider for this? edit: To be clear, I am not looking at pairwise relationships. What I am looking for is a method that would look at the mountain of data at hand and attempt to discover (potentially non-linear) lead-lag relationships between _arbitrary_ groups of series and the individual series to be predicted.


title: Confidence Bands in non linear mixed effects models

text: Is there a way to get confidence bands around the fitted population fixed effects curve for non linear mixed effects models fitted in `R` function `nlme()`? I would like to plot the population mean fitted curve with the uncertainty in the estimates, however, I have not found any clear way to do this taking into account the random effects variance.


title: Sample size for multiple regression: How much more data do I need?

text: I'm trying to test a prediction from a computational model with a biological experiment. The model predicts a linear dependence between cell size and asymmetry. I performed an experiment that is known to affect cell size, and got a nice significant slope when regressing asymmetry onto cell size. However, in biology, there is rarely a perfectly clean experiment. To account for potential other effects of the perturbation, I added a dummy variable for the perturbation to the regression. The coefficient of the dummy is not significant (p=0.5), but the coefficient of cell size is no longer significant, either (p=0.1). Assuming that the model prediction is correct, how can I, given the results of my regression, calculate how many additional measurements I need to perform so that the coefficient of cell size remains significant even in the presence of the confounding dummy?


title: Confusion regarding correlation and mse

text: I have a dataset containing three variables. I have this confusion, low mse means higher correlation isn't it? Now when I calculate the mse for variable two and three where the third variable is my ground truth, I get higher mse than the one I get by comparing variable one and variable three (ground truth). However, when I calculate the correlation using corr function in matlab, I get higher correlation between variable two and three around .78. On the other hand, the correlation between variable one and three is .76. Why is it so? Let me tell you I have different number of variable1-variable3 and variable2-variable3 set. I meant to say sometimes the values of the variables are missing. So altogether, I have 1280 sets of variable1-variable 3 pairs and 400 of variable2-variable 3 pairs. As I have said I have higher mse for the variable1-variable 3 set than variable2-variable3 set calculating for the number of pair points 400 and 1280 respectively. I am a bit confused. Any insights?


title: Simple linear regression fit manually via matrix equations does not match lm() output

text: I am trying to fit a linear model using matrices to my data set even though I can use OLS and do it without matrices as a simple tutorial for myself to better understand both `R` and matrix notation. This is the model I am trying to fit: $$\bf Y=X\boldsymbol\beta+\varepsilon$$ where $\bf Y$ is a $1\times n$ matrix, $\bf X$ is a $n\times k$ matrix (where $k$ is the number of $\beta$'s, which in this case is 2), $\boldsymbol \beta$ is a $k\times 1$ matrix and lastly our error term is $n\times 1$. I understand this portion. When I simply use the `lm()` command to fit my data, I get the following from the `summary()` command:               Call:     lm(formula = y ~ x)          Residuals:         Min      1Q  Median      3Q     Max      -4.0503 -1.4390  0.4921  1.0589  3.9446           Coefficients:                 Estimate Std. Error t value Pr(>|t|)         (Intercept)   2.9849     0.8219   3.632  0.00191 **      x             0.5612     0.1084   5.178 6.32e-05 ***      So the `summary()` is telling me that the $\beta$ matrix is a $2\times 1$ matrix, with the first number (which is $\beta_0$) as 2.0949 and the second number (which is $\beta_1)$ as 0.1084. My question is this: We know that the matrix $\beta$ is actually: $$\boldsymbol\beta=(\bf{X}^T\bf{X})^{-1}(\bf{X}^T\bf{Y})$$ and when I try to simply carry out this calculate by hand using R using `b=(t(x)*x)^-1*t(x)*y`, I get a $1\times 20$ vector (where $20$ of course is $n$, the number of observations). Why am I not getting a $2\times 1$ matrix like I should be getting?


title: Convergence of MCMC for ill-behaved functions

text: We are doing MCMC sampling over an n-dimensional space with a population of MCMCs. The goal is to have them stop when they approximately converge. The problem is that the convergence, which is so easily visually seen is hard to code. Typical Gelman-Rubin/Geweke work well for well behaved functions, but for some nasty ones with weird features, they oscillate or have random peaks that do not indicate the steady slow convergence as wanted. What other convergence statistics can you recommend? I was looking at using time series methods but cannot find any papers on time series convergence. I was thinking would it be possible to take a discrete fourier and do something with that or.....? Some already looked at: Cramer von Mises, Kolmogorov-Smirnov, Kuiper, Raftley- Lewis, Heinberg-Wurtley.


title: Whether to use original or reverse coded items in factor analysis?

text: I am currently analyzing data from a 34-item Likert scale. I already recoded the negatively stated variables in SPSS as different variables. I'm about to do a factor analysis. **Should I use the original variables or the recoded variables in the factor analysis?**


title: Advice on sensitivity analysis for priors in Bayesian statistics

text: I'm not clear on how to perform sensitivity analysis on the priors. Many sites have different answers. One site indicates to perform three non-informative, weakly informative and known priors. Another suggest running the model with different priors. Here are my questions: I want to estimate this parameter A, in which A = mean1 + mean2 . So, the priors for mean1 and mean2 are formulated ( **in coding** ) as the following :               mean1 ~ N(u1, tau1)     A ~ N(mean1, tau2)       where               tau1 ~ gamma(alpha1,beta1)     tau2 ~ gamma(alpha2,beta2)      With the above formulation, how do I perform the sensitivity analysis?   1. Should I define a range of u1, tau1, and tau2 ? If so, what is the systematic procedures? It seems to me there are various combinations: varying u1 and fixing tau1 and tau2; or fix u1 and tau2, varying tau1; or fix u1 and tau1, varying tau2. How do we know what are the correct priors to use?    2. How do I incorporate weakly informative priors to the above formulation?    3. What if the sensitivity test returns different estimated mean1 and A by varying the priors? How do we make conclusion? Your advice will be appreciated and help me get started. Thank you! EDIT: I corrected an error. You can calculate the estimated _mean2_ from estimated parameters _A_ and estimated _mean1_.


title: similarity measures with missing values

text: How to handle missing values when computing similarity (or distances)? (I have binary feature values and do use the simple matching coefficient, but I feel that the answer to this question may be more general) I can think of two options:   * Remove missing values   * Count missing values as error But removing null values has the problem that a high score can be achieved with only one/a few values (see example A). And counting missing as error has the problem that real miss matches should be counted higher than missing values (see example B). Is there a technique that has none of these shortcomings? Example A                          Instance1 Instance2       Feature1 missing   missing       Feature2 1         missing       Feature3 0         0       Feature4 missing   1       Feature5 missing   missing       Feature6 missing   missing            Simple-matching-similarity-REMOVE = 1/1 (twice as high as B)       Simple-matching-similarity-COUNT-AS-ERROR = 1/6      Example B                          Instance1 Instance2       Feature1 missing   1       Feature2 1         0       Feature3 0         0       Feature4 1         1       Feature5 1         0       Feature6 0         missing            Simple-matching-similarity-REMOVE = 2/4       Simple-matching-similarity-COUNT-AS-ERROR = 2/6 (twice as high as A)


title: Summing variance of autocorrelated timeseries

text: **Problem** : When trying to calculate the variance of timeseries sums I get a negative variance, mostly due to autocovariances at large lag steps. Does not seem realistic. I have a timeseries which is calculated from another timeseries using a regression equation. I would like to propagate the uncertainty in the regression to the final timeseries. Then I want to sum (or take mean values) different segments of the timeseries over different timeperiods, and get the uncertainty of the sums. The timeseries is originally in 1 hour frequency and I want to sum over periods of 1 day (resampling to daily frequency) up to several years. The timeseries is strongly autocorrelated at short lag times. For getting the variance of the sum (in the case of 3 elements being summed): $$Var(a+b+c)= \\\ Var(a)+Var(b)+Var(c) + 2 \times (Cov(a,b) + Cov(a,c)+Cov(b,c))$$ I use `r` for the calculations. I get the variances for each timeseries element as $SE^2$, where $SE$ is the standard error (`se.fit`) returned from r's `predict()` function using the regression model. The covariances I get from the autocovariance function `acf()`. Here is some code and a selection of the data (excuse clumsy R code, I'm very new to R):               #tsY is the predicted timeseries from the regression     tsY=c(81.4,  79.0,  83.4,   81.7,   75.7,   68.3,   62.3,   57.2,   52.6,   48.8,   45.4,   42.6,   39.9,   37.6,   35.6,   33.8,   32.2,   30.8,   29.6,   28.4,   27.3,   26.2,   25.0,   23.9)     #tsSE is the standard error from the prediction (se.fit)     tsSE=c(1.55,  1.49, 1.60,   1.56,   1.41,   1.23,   1.09,   0.97,   0.87,   0.78,   0.71,   0.65,   0.60,   0.55,   0.51,   0.48,   0.45,   0.42,   0.40,   0.38,   0.36,   0.34,   0.32,   0.30)          tsVar=tsSE^2          #create a matrix of the autocovariances at different lag times, diagonal is lag=0     #rows and columns are indicies in timeseries     covmat<-matrix(numeric(0), length(tsY),length(tsY))      for ( i in (1:(length(tsY)) ) ) {       if (i == 1) {         autocov<-acf(tsY, type='covariance', lag.max= length(tsY))         autocovvec<-autocov$acf[1:nrow(autocov$acf)]         covmat[i:length(tsY),i]=autocovvec       }  else {         autocov<-acf(tsY[-(1:i-1)], type='covariance', lag.max= length(tsY))         autocovvec<-autocov$acf[1:nrow(autocov$acf)]         covmat[i:length(tsY),i]=autocovvec       }          }          # sum the matrix columns, but not the diagonal     sumofColumns <- rep(NA, ncol(covmat))     for (i in (1:ncol(covmat))) {       if (i == 1) {         sumofColumns[i]=sum(covmat[-(1),i])         } else{          sumofColumns[i]=sum(covmat[-(1:i),i])         }     }          sumofCov=sum(sumofColumns) # sum of the covariance (Cov(a,b) + Cov(a,c)+...)     sumofVar=sum(tsVar) # sum of the variances of each timeseries element     varofSum=sumofVar+2*sumofCov # variance of the sum of the timeseries          # from the covmat the negative variance occurs at larger lag times.     acf(tsY, type='covariance', lag.max= length(tsY))          > sumofCov     [1] -1151.529     > varofSum     [1] -2283.246      **So I have the following questions:** >   1. Did I completely misunderstand how to calculate variance of sums? > >   2. Is it better to use a cutoff from the max lags to be considered in the > autocovariance? If so how would one determine this? This would especially be > important with the complete data where the length is several thousand. > > > > ~~3\. Why is the covariance negative in this sample data at large? When > plotting tsY`plot(tsY)` it looks like the covariance/correlation should > remain positive.~~ Because it is the variation in direction from their > means. **EDIT:** > Comment on **question 2** above: I have realized that using n-1 lags, as > above in the code, does not make a lot of sense. There appear to be few > different ways to determine the maximum lags to consider. Box & Jenkins > (1970) suggest n/4 and R by default 10*log10(n). This does not answer the > question however, of how to determine an appropriate cutoff for summing the > covariances. > > Does it make sense to look at the partial autocorrelation (function pacf()), > in order not to overestimate the effect of the auto covariance in the > summation term? The partial autocorrelation for my data is significantly > different from zero only at 1 or 2 lags. Similarly, fitting an AR model > using ar() function, I also get an order of 1 or 2. Cheers Related post Variance on the sum of predicted values from a mixed effect model on a timeseries


title: R modeest package: Which method of mode estimation?

text: I have a unimodal continuous distribution which is skewed and I would like to estimate the mode using the modeest package in R. However, the package offers a relatively large number of methods to do so. What makes the situation more complicated is that some of these methods have even more options for different ways of estimating the mode within the particular method. Since I'm a relative newbie to statistics in general, I'm a little overwhelmed. How can I choose the right method here? Is it a valid approach to just plot the distribution and use the method which returns a value that seems reasonable given the graphic information?