title: Model selection with AIC and model efficiency

text: Starting from an "a priori" set of models based on my knowledge about potential relations between my dependent variable and the independent variables considered, I use AIC to find best models. More specifically I calculate Akaike weights then Evidence Ratio (ER) and consider that models with a ER < 2 are equally likely. But the same problem remain each time I do that. I selected the best models from a set of them, but I don't know if those models are efficient to predict (or at least represent) my data. I can have selected the best element(s) of the list of the worst models. I do not use $R^2$ in model selection because of the fact that including more variables generally increase $R^2$ value. But ! When the selection by AIC is done and I can consider that models with a ER < 2 are equally likely. Do you find it is correct to calculate $R^2$ or pseudo-$R^2$ for the best "set of models" in order to have an idea of the representativeness of those models and use this value to select the more efficient model? I would be glad to hear your opinions about this! Note: Thanks for suggestions about cross-validation, I will try this. Unfortunatly, I do not have an external dataset to test my models.


title: Ripper algorithm on large data sets

text: I've got a supervised learning problem where I've got about 15 features, mostly numeric, and I'm mapping these to a set of 5-7classes. Using decision trees or forests, I'm able to get confusion matrices that I'm quite satisfied with, but I'm interested in seeing what the "top" rules RIPPER would find would be. I've been attempting to train on a small subset of the data, about 2 million rows (Again 15 features), but I find myself unable to actually generate a model via JRip as I always run into memory constraints after 6-7 hours of Weka's JRip churning. (Even on large EC2 instances) Has anyone had success using RIPPER on large-ish data sets? Would I find more success 'bucketing' my numeric attributes first? Are there any important preprocessing steps I'm missing?


title: How do I interpret interaction coefficients in a temporally lagged regression model?

text: I am running a linear regression on daily price returns and I want to capture the interaction between subsequent returns: $r_{t+1} = \alpha + \beta_1 r_t + \beta_2 r_{t-1} + \beta_3 r_t \cdot r_{t-1} + \epsilon$ Where $r$ is a geometric return ($p / p_{t-1}$). Each coefficient is interperable as a linear multiplicative factor of the return. How can I interpret the meaning of $\beta_3$, the interaction term?


title: Optimal strategy in growing jackpot scenario

text: There is a betting game, which seems to have "positive sum" situations, and I'd like to find out the optimal strategy. Multiple players place bets against the bank, for and bets that fit winning criteria earn players points of the game. There are two sets of criteria: The first: bet 1 point with 0.475 chance of doubling it or 0.525 – losing it. This is obviously a negative sum game. The second part is the "Jackpot": you can add a fixed amount: 2pt, to your prior bet, to have the bet included in roll for the jackpot. The chance for winning the jackpot is 0.0001. It's all about how big the jackpot is. All points you lose, be it to the 0.525 chance of loss in the first part, or the 0.999 chance of loss in the second part, are accumulating in the jackpot. Additionally, if the amount in the jackpot is less than 10,000, the house covers the difference between the amount and 10,000 (you can't win less than 10,000 if you hit the 0.0001 chance). When a winning bet for the jackpot is placed, the jackpot is emptied, and the process of collecting begins anew. When always playing the first part as bet 1, with house edge of 5%, and jackpot bet cost of 2, the expected value for a bet is (-2.05 + (jackpot × 0.0001)) which means it becomes a positive sum game when jackpot value exceeds 20,500. When that happens I should keep placing bets until I win the jackpot. But I don't quite grasp my chances for jackpot in the two lower ranges: 0–10,000 and 10,000–20,500. If it was a simple case of a fixed-amount jackpot, that would be plainly a negative sum game, but considering all my bets increase the jackpot value, I'm not so sure. Also, does the fact that the jackpot is shared between all players (their losses add to it, one takes it all) change the optimal strategy?


title: Why does chi-square testing use the expected count as the standard deviation?

text: In $\chi^2$ testing, what's the basis for using the square root of the expected counts as the standard deviations (i.e. the expected counts as the variances) of each of the normal distributions? The only thing I could find discussing this at all is http://www.physics.csbsju.edu/stats/chi-square.html, and it just mentions Poisson distributions. As a simple illustration of my confusion, what if we were testing whether two processes are significantly different, one that generates 500 As and 500 Bs with very small variance, and the other that generates 550 As and 450 Bs with very small variance (rarely generating 551 As and 449 Bs)? Isn't the variance here clearly not simply the expected value? (I'm not a statistician, so really looking for an answer that's accessible to the non-specialist.)


title: Designing a Simple (a.k.a 'bad') Ranking value from several values of unknown distribution

text: I am but a humble developer, so please forgive if I butcher terminology of confuse concepts. I am looking at designing a simple ranking of **Average Internet Audience Rank** based upon a bloggers audience size (Twitter/Tumblr/Facebook friends) and how many times their content has been shared (i.e. 'Reach'). **Here are the metrics I want to combine:**               | metric       | weight |    range      | distribution (est.)     -------------------------------------------------------------     | shares       |  1st   |  0 - Millions | power law     | comments     |  2nd   |  0 -thousands | power law     | Klout score  |  3rd   |  0 - 100      | normal , weighted toward lower end     | Twitter fol. |  4th   |  0 - Millions | power law     | LinkedIn con.|  5th   |  0 - 1000     | normal      **Goal**   * None of the individual metrics will completely crowd out any of the others if it is significantly higher than normal. For example if someone had 10,000,000 twitter followers, but all other values were 0, I would never want them to rank higher than someone who had even 1000 blog post shares. I know this Multi-Criteria Decision Analysis which I've been reading up on, but operations research is not a field I'm very familiar with. Hopefully you guys can help point me in the right direction.


title: Edge correction of Ripley's K-function for two 1D point processes?

text: I am just beginng an investigation involving characterizing the dependence between two 1D stochastic point processes $x$, $y$. The natural approach seems to involve Ripley's K-function: $$ K(t) = \frac{T}{n_xn_y} \sum_{i=1}^{n_x} \sum_{j=1}^{n_y} w(x_i,y_j) I[d(x_i,y_j)<t] $$ where $n_x$ is the number of observations in $x$ and $n_y$ is the number of observations in $y$ across the interval $T$. Deviation from $K(t)=t $ is an indication of correlation between the two point processes. However, it's not clear to me how to estimate the edge correction $w(x_i,y_j)$ the 1D case. Some papers refer to Hani Doss' 1989 papers, but in the JSTOR paper he explicity states that ' this edge correction will not concern us.' FWIW-I am currently using a correction weight of 2, but intuitively that seems excessive in my case.


title: Handle leap seconds in time series data

text: I am trying to figure out how to handle time series data. Initially, I thought it was self-evident (for my purposes) to only allow non- overlapping data blocks. (I.e., non-overlapping, time wise.) Then it occurred to me, what about leap seconds? There could be a leap second in there, making one of the hours 3601 seconds long instead of 3600, and throw the whole thing off. How do people handle things like these? Shaving off the last second? Have a "fuzz" factor and allow a second up or down? Am I over-thinking this?


title: Safely determining sample size for A/B testing

text: I am a software engineer looking to build an A/B testing tool. I don't have a solid stats background but have been doing quite a bit of reading over the last few days. I am following the methodology described here and will summarize the relevant points, below. The tool will allow designers and domain experts to configure a website to split traffic received at a specific URL between two or more URLs. For example, traffic arriving at http://example.com/hello1 could be split between http://example.com/hello1 and http://example.com/hello2. Traffic would be split evenly between target URLs and the performance of the marketing processes at each of the target URLs will be compared. In this experiment, the sample size `N` will correspond to visitors. The test will measure "conversions", a term describing when a visitor commits to a specific action in a marketing process. Conversions are expressed in percentages and a higher conversion rate is desirable. This makes the test a comparison of independent proportions. The tool needs to be able to be easily employed to produce tests with safe results. Selecting an appropriate value of `N` is important. In the linked article, above, a power analysis of two independent proportions is employed to find `N`. This method requires that one know the conversion rate of the control in advance as well as specify the target desired conversion improvement. It also specifies a significance level of 95% and a statistical power of 80%. Questions:   1. Is this method of determining `N` sound? If so, what is the safest way to determine the conversion rate of the control prior to beginning the test?   2. Are there sound ways of determining `N` that don't require that one know conversion rates of the control in advance?   3. Is the methodology in the linked article sound? If not, are there any accessible and easily digestible methods out there that you could link me to?


title: Standard error from OR and HR

text: I'm planning to combine several observational studies, some of which have produced odds ratios and some hazard ratios. How do I get the standard error from the odds ratio $1.57$ ($95\%$ CI: $1.08-2.27$) and the standard error from the hazard ratio $1.56$ ($95\%$ CI: $0.96-2.52$)?


title: sample size for well-established questionnaire

text: I would like to conduct an usability study to assess users' perceived satisfaction with a computer system. I know that there are many well- established questionnaires to carry out this assessment. Now I am wondering how I can determine the sample size to be used in my study.Is there a specific rule?


title: Resource to read about distributions

text: My course textbook for statistics is the one by Levine and Smith. However, I find myself interested in the various probability distributions. The author simply states the $f_X(x)$ and derives $F_X(x)$ and gives an application. I was in the hunt of a good online ( _FREE_ ) resource which talks about probability distributions in a more intuitive setup. For instance, telling me WHY the equation for exponential, Weibull, Normal, Gamma etc. are the way they are. Especially gamma distribution, he simply states that in the calculus class one can learn to prove that the integral is convergent and states $f_X(x)$.


title: Radial profile and 2d log-normal distribution?

text: I have a case study where a person should be located. We do not know where this person is, but we have some information. The total story is basically about the person which is to be searched. The informations we have are connected to a road, a house and a satelitte: One part is: > A probability distribution centered around a house also informs us of the > person's location. The distribution's radial profile is log-normal with a > mean of 4800m and a mode of 3900m in every direction. My question: **How to model this log-normal distribution?** If I am correct, the information about the radial profile tells me that I should model it with a 2d log-normal distribution? Is this right? How is the 2d probabiliy distribution of a log normal distr. and how do I model it in R?


title: Predict best action sequence for dataset of debt observations

text: I need some suggestions of general predicting practice in the following case: I have a dataset of debt observations; there are ~8 variables defining each debtor situation (debt details and person details), quantitative and qualitative. The Debt Management applied for each person **a sequence of actions** (there are a few types of actions, each of different cost and different average effectiveness). The results I have in use consist of 2-elem vector:               [agreement, success_in_court]        * `agreement` is of `1/0` value, stands for reaching an agreement with debtor or not   * `success_in_court` is of `1/0` value, stands for reaching success in court in case this action (to sue a debtor) was applied by the Debt Management  Logically, no more than one element of a "response vector" may be of value 1. The point is to predict the most profitable sequence of action for future debtors. What methods should I use? It sould be implemented in `R`.


title: Given a sample of time-series spatial data, how do I determine the minimum number of spatial points that is statistically similar to the large set?

text: I have a solar irradience data set that was generated from 45 independent stations within a square mile. The stations are more or less regularly spaced in a grid. The data are 1Hz and are all time-synced (in theory). Having all 45 stations gives us a good idea of what the sunlight is doing over this square mile for each second of the 550 days of data that we have (20Gb). However, photovoltaic developers don't want to spend money on 45 stations to measure solar irradience when they can spend money on X stations over Y area and determine the solar variability. The question is what are X and Y? Since maintaining the solar variability is necessary, the area is part of the equation. Solar variability refers to changes in measured solar irradience over space and time. Clear and overcast days are not variable whereas partly cloudy days can generate larges spikes in irradience values. The spikes in irradience throughout the day are quantified as having high "solar variability". For PV power plants the spikes result in voltage changes on transmission lines, something utilities do not like. For more information see: http://energy.sandia.gov/wp/wp- content/gallery/uploads/Stein_ASES_2012_VI_paper_SAND2012-2088C1.pdf I wrote a Matlab program to tackle this problem. First, I used 10 seconds of the full data set and calculated the mean and StdDev for each second using all 45 stations. (Once the algorithm works on these 10 seconds, I'll use the larger set.) Then I randomly chose a sample subset (with replacement) of N elements from the 45 stations for each second where 3 < N <=15. I chose the N elements randomly about 100 times each and calculated the average of all the 100 standard deviations. Then I looked for an N that had an average standard deviation value close to "the" answer for the 45 stations. As it turned out, when N=14 then the StdDev of the subsets closely matched "the" set. I also looked at the area that those random samples enclosed (the convex hull) though results have been inconclusive. I calculated metrics such as 1) area enclosed 2) perimeter of convex hull 3) sensor density for area 4) sensor density for perimeter. I'd just like to know if I'm on the right track here, or if there's a method that I don't know about that would be a good fit for this application.


title: Conditional expectation of R-squared

text: Consider the simple linear model: $$\pmb{y}=X'\pmb{\beta}+\epsilon$$ where $\epsilon_i\sim\mathrm{i.i.d.}\;\mathcal{N}(0,\sigma^2)$ and $X\in\mathbb{R}^{n\times p}$, $p\geq2$ and $X$ contains a column of constants. My question is, given $\mathrm{E}(X'X)$, $\beta$ and $\sigma$, is there a formula for a non trivial upper bound on $\mathrm{E}(R^2)$*? (assuming the model was estimated by OLS). *I assumed, writing this, that getting $E(R^2)$ itself would not be possible.  # EDIT1 using the solution derived by Stéphane Laurent (see below) we can get a non trivial upper bound on $E(R^2)$. Some numerical simulations (below) show that this bound is actually pretty tight. Stéphane Laurent derived the following: $R^2\sim\mathrm{B}(p-1,n-p,\lambda)$ where $\mathrm{B}(p-1,n-p,\lambda)$ is a non-central Beta distribution with non-centrality parameter $\lambda$ with $$\lambda=\frac{||X'\beta-\mathrm{E}(X)'\beta1_n||^2}{\sigma^2}$$ So $$\mathrm{E}(R^2)=\mathrm{E}\left(\frac{\chi^2_{p-1}(\lambda)}{\chi^2_{p-1}(\lambda)+\chi^2_{n-p}}\right)\geq\frac{\mathrm{E}\left(\chi^2_{p-1}(\lambda)\right)}{\mathrm{E}\left(\chi^2_{p-1}(\lambda)\right)+\mathrm{E}\left(\chi^2_{n-p}\right)}$$ where $\chi^2_{k}(\lambda)$ is a non-central $\chi^2$ with parameter $\lambda$ and $k$ degrees of freedom. So a non-trivial upper bound for $\mathrm{E}(R^2)$ is $$\frac{\lambda+p-1}{\lambda+n-1}$$ it is _very_ tight (much tighter than what I had expected would be possible): for example, using:               rho<-0.75     p<-10     n<-25*p     Su<-matrix(rho,p-1,p-1)     diag(Su)<-1     su<-1     set.seed(123)     bet<-runif(p)      the mean of the $R^2$ over 1000 simulations is `0.960819`. The theoretical upper bound above gives `0.9609081`. The bound seems to be equally precise across many values of $R^2$. Truly astounding! # EDIT2: after further research, it appears that the quality of the upper bound approximation to $E(R^2)$ will get better as $\lambda+p$ increases (and all else equal, $\lambda$ increases with $n$).


title: Weighted average of measurements with unequal errors

text: Suppose I have the numbers below. Consider them as results of some measurement.                22.73  +2.15  -1.81      23.81  +1.77  -1.54      28.6  +11.4   -6.4       29.59  +8.58  -5.43      The numbers on the 2nd and 3rd columns are the errors associated to the values in the 1st one. I would like to know how to calculate an average of those numbers taking the errors in account. According to what I saw on the net, I should use a weighted average, but I was unable to find an example which uses values like those above, i.e. with unequal errors.


title: What is the variance of the mean of correlated binomial variables?

text: An average of $B$ binomial i.i.d. random variables, each with variance $\sigma^2,$ has variance $\frac{1}{B}\sigma^2.$ If the variables are simply i.d. (identically distributed, but not necessarily independent) with positive pairwise correlation $\rho$, the variance of the average is $$\rho\sigma^2 + \frac{1-\rho}{B}\sigma^2$$ but I don't understand why. Can somebody provide a proof?


title: Prognostic factor vs. risk factor

text: As far i understand the difference between `Risk factor` and `Prognostic factor` is : Risk factor is the `cause` of a condition whereas Prognostic factor influences the outcome who has the condition, like resectability of tumour for lung cancer. Is that the case ? Would you explain the difference with a similar example in both cases ?


title: Finding a correspondence between time-series elements

text: My problem deals in particular with time-series data about server performance, but the solution is sure to be applicable to many types of data sets. Pardon me if the answer is well-known; I don't know the right vocabulary to effectively search for these things. Suppose I have two sets of time-series observations (let's call them A and B) containing thousands to millions of observations each. For example,               A: [ (1367678860, 30), (1367678870, 40), (1367678880, 33), ... ]     B: [ (1367678863, 0.18), (1367678873, 0.11), (1367678884, 0.12), ... ]      The data in `A` and `B` are, in general, collected by different processes with different sampling rates. Furthermore, any number of data points could be straight-up missing from either series. I want to _parameterize_ this data by time with a one-to-one correspondence. That is, I need an algorithm that will find a sequence of ordered pairs `(a_i, b_i)` whose corresponding times are close to each other. For example, the data above would be converted to:               [ (30, 0.18), (40, 0.11), (33, 0.12), ... ]      As I said, it needs to work on data sets that don't necessarily share a sampling rate, and for which arbitrary points may be missing. Points that don't match up nicely (whatever that means) with the other time series should be discarded. I've been banging my head against this for a while, trying to adapt the least- squares linear regression method to the problem, but so far I haven't come up with anything good. Stats nerds: what do you suggest?


title: Within-subject and between-subject fixed effects in mixed model

text: I've been trying to analyze some data using mixed models but I have some troubles to understand how should I include both within-subject and between- subject fixed effects in such models. Let's consider a following theoretical example in which I have pupils nested within schools. I want to analyze the effect of the class size on the pupil's performance (in 1-10 scale), BUT also the effect of the fixed effects not nested within each subject (school), such as the school size:               school   school_size   pupil   class_size   performance     S1       500           1       30           5     S1       500           2       20           6     S1       500           3       20           8     ...     S20      1000          101     20           4     S20      1000          102     30           3     S20      1000          103     30           10      Is the following mixed-model appropriate in this case?               performance ~ school_size + class_size + (1|school)      I particularly would like to know if I can include both within-school (i.e. class size) and between-school (school_size) fixed effects, or should I analyze between-school effects with linear model and within-school effects with a second, mixed model?


title: Weighted principal components analysis

text: After some searching, I find very little on the incorporation of observation weights/measurement errors into principal components analysis. What I do find tends to rely on iterative approaches to include weightings (e.g., here). My question is why is this approach necessary? Why can't we use the eigenvectors of the weighted covariance matrix?


title: Assigning class labels to k-means clusters

text: I have a very basic question on clustering. After I have found k clusters with their centroids, how do I go about interpreting the classes of the data points that I have clustered (assigning meaningful class labels to each cluster). I am not talking about validation of the clusters found. Can it be done given a small labelled set of data points, compute to which cluster these labelled points belong to and based on type and number of points each cluster receives, decide the label? This seems pretty obvious but I don't know how standard it is to assigns labels to clusters this way. To be clear, I want to perform unsupervised clustering that doesn't use any labels to first find my clusters. Then having found the clusters, I want to assign meaningful class labels to the clusters based on the properties of a few example datapoints.


title: log-likelihood and negative binomial zero inflated distribution in glmmADMB

text: Anyone know how glmmADMB calculates log-likelihood and negative binomial zero inflated distribution? I compared a model using glmmADMB and pscl and the diference in AIC and parameters are huge. I don´t know why. HELP.


title: Is the exact value of a 'p-value' meaningless?

text: I had a discussion with a statistician back in 2009 where he stated that the exact value of a p-value is irrelevant: the only thing that is important is whether it is significant or not. I.e. one result cannot be more significant than another; your samples for example, either come from the same population or don't. I have some qualms with this, but I can perhaps understand the ideology:   1. The 5% threshold is arbitrary, i.e. that p = 0.051 is not significant and that p = 0.049 is, shouldn't really change the conclusion of your observation or experiment, despite one result being significant and the other not significant.  The reason I bring this up now is that I'm studying for an MSc in Bioinformatics, and after talking to people in the field, there seems to be a determined drive to get an exact p-value for every set of statistics they do. For instance, if they 'achieve' a p-value of p < 1.9×10-12, they want to demonstrate HOW significant their result is, and that this result is SUPER informative. This issue exemplified with questions such as: Why can't I get a p-value smaller than 2.2e-16?, whereby they want to record a value that indicates that by chance alone this would be MUCH less than 1 in a trillion. But I see little difference in demonstrating that this result would occur less than 1 in a trillion as opposed to 1 in a billion.   2. I can appreciate then that p < 0.01 shows that there is less than 1% chance that this would occur, whereas p < 0.001 indicates that a result like this is even more unlikely than the aforementioned p-value, but should your conclusions drawn be completely different? After all they are both significant p-values. The only way I can conceive of wanting to record the exact p-value is during a Bonferroni correction whereby the threshold changes due to the number of comparisons made, thus decreasing the type I error. But even still, why would you want to show a p-value that is 12 orders of magnitude smaller than your threshold significance?   3. And isn't applying the Bonferroni correction in itself slightly arbitrary too? In the sense that initially the correction is seen as very conservative, and therefore there are other corrections that one can choose to access the significance level that the observer could use for their multiple comparisons. But because of this, isn't the point at which something becomes significant essentially variable depending upon what statistics the researcher wants to use. Should statistics be so open to interpretation?  In conclusion, shouldn't statistics be less subjective (although I guess the need for it to be subjective is as a consequence of a multivariate system), but ultimately I want some clarification: can something be more significant than something else? And will p < 0.001 suffice in respect to trying to record the exact p-value?


title: Elementary MCMC pseudocode

text: My aim relates to the project described here, but I've tried to make this question self-contained. I'm trying to write the MCMC pseudocode for the following inference problem: > Given two observed independent outcomes $X$ and $Y$, what are the posterior > densities of probabilities $p_{x}$ and $p_{y}$ of obtaining $X$ and $Y$, > assuming $X$ ~ $\text {Bin}(N, p_{x})$ and $Y$ ~ $\text{Bin}(N, p_{y})$? We > do not know or observe $N$. We have priors for $N$, $p_{x}$, and $p_{y}$. We > will assume that $N$ ~ $\text{Poisson}(\lambda)$, and $p_{x}$ and $p_{y}$ > will have either $\text{Beta}$ or $\text{Uniform}$ priors, chosen by the > user. For the $\text{Beta}$ priors, the mean and variance will be specified. > For the $\text{Uniform}$ priors, the upper and lower limits will be > specified. My background is obviously not statistics, and I'm having a hard time patching gaps in my knowledge from the voluminous MCMC literature. Rather than burden commenters with explaining a mountain of Bayesian inference to me, **I'd appreciate it if people could point out the errors in my pseudocode and understanding, including instances in which I might be misusing terms**. I'll focus my study on those areas. I welcome suggested reading that is specific (I've seen but not completely read some of the suggestions here). I will be coding this by hand from scratch, except for simple math functions. Ultimately, I'll make sure there's complete and accurate pseudocode associated with this post. **Proposed pseudocode** >   1. Declare an initial value of $\theta_{0} = (N,p_{x},p_{y})$ and set > equal to $\theta_{\text{current}}$. >   2. Calculate the posterior probability of $\theta_{\text{current}}$: > >     * First, calculate the log-likelihood $\textit{L}_\text{current}$ of > $\theta_{\text{current}}$. Here, the log-likelihood > $\textit{L}_\text{current}$ = $L_{x} + L_{y}$, where $L_{x}$ $=$ > $\text{log}{N \choose X} + X$ $\text{log } p_{x} + (N-X) \text{ log }$ > $(1-p_{x})$ and $L_{y}$ is analogous. > >     * Next, obtain the posterior by multiplying the log-likelihood > $\textit{L}_\text{current}$ by the prior, > $\text{Prior}(\theta_\text{current})$. The prior of $\theta_\text{current}$ > is the product of the priors of each of the parameters in $\theta$. In this > case, if the user has chosen $\text{Beta}$ priors for $p_{x}$ and $p_{y}$, > the respective priors are simply the $\text{Beta}$ pdf evaluated at $p_{x}$ > or $p_{y}$. The prior for $N$ is the pmf of the Poisson function at $N$. > >   3. Draw $\theta_{\text{candidate}}$ from the proposal distribution. **What > exactly is this distribution for each parameter in $\theta$?** I've seen > normal distributions used, but here $N$ and $p_{x}$ clearly must be drawn > from different distributions. >   4. Calculate the posterior of $\theta_{\text{candidate}}$ following step > 2. >   5. Evaluate whether to accept $\theta_{\text{candidate}}$ by comparing the > posteriors, and update $\theta_{\text{current}}$ if justified. I will use > the Metropolis algorithm (not shown for brevity). >   6. Repeat steps 3-5 until convergence or reasonable stopping point. > **(tl;dr) Major questions**   * In step 2, am I missing shortcuts related to the conjugate priors? I still don't really understand what they're good for or if I'm calculating the prior correctly.   * How are jump distributions determined? Thanks in advance.


title: Missing values in GLM

text: I am applying glm on a data in which most of the values are NAs or blank. For example, in the example data produced below (4 predictors and one response variable), the default glm command will remove 10 rows that contain 'NA' leaving just one row for analysis. This creates serious problem as some of my data that initially had 100000 rows (with ~50 features) cut down to 200 with (~15 features) or even less reducing the power significantly. My question is: What options do I have in this scenario. I do not want to fill the NAs with average values or anything based on variance/measure of centrality as it might be possible that TSH test is ordered only for patients with history of thyroid disease. In that case any extrapolation will be disastrous, as the mean won't be the actual representatives of the normal values.               gender  TSH PH  HDLC_hole   response     m   2   NA  36  TRUE     f   1.8 4   32  TRUE     m   NA  NA  29  TRUE     f       NA  NA  33  TRUE     m   2.2 5   NA  TRUE     f   2.5 4   NA  TRUE     NA  1.8 4   34  FALSE     m   NA  4   35  FALSE     f   3   NA  36  FALSE     m   1.2 4   NA  FALSE     m   1   NA  28  TRUE      Edit: What I understood from a quick review of some of the papers on 'Multiple imputation' techniques is that they fill-in values based on mean/variance or similar statistics. Should I change my method and explore others? For such cases what predictions algorithms may be appropriate.


title: Tricks for a very fast implementation of Random Forest

text: I am implementing my own Random (regression) Forest algorithm and am looking for tricks to speed up the estimation of forests on large datasets. So far I have implemented three main tricks: 1) Use a running variance on the sorted response vector to efficiently determine the optimal split point (for regression forests). 2) Determine whether a split is optimal only on a subset of all possible split points. 3) Parallelize the construction. What other tricks are there that I am missing out on? Note that CPUs are scarce but hard drives are not, so I am looking for tricks like some kind of pre-sorting or other such tricks that save information to disk that can be loaded whenever I estimate a new tree to speed up that tree's estimation. RAM is somewhat scarce but it's mostly CPU scarcity since I can use multithreading over a large RAM base/heap size. I am sure that my implementation can be sped up by 3-4x but I am not sure how.


title: Paired t-test and two-sample t-test

text: I have a simple question regarding the use of paired t-test and two-sample t-test:   1. What is the consequence of using the independent two-sample t-test for dependent paired samples?   2. What is the consequence of using the paired t-test for independent samples? Thanks!


title: Non-response bias (many waves, time span between waves not consistent)

text: I want to ask regarding non-response bias. I have sent my questionnaire on- line to SMEs in Malaysia. However, I have sent reminders many times (around 5 times). However, now I cannot recall when the reminders were taken place as my account has been cancelled on the on-line questionnaire system. However, based on my guessing, I am doing ONE way ANOVA to see if there is a significant difference between the data collected based on the waves. the test shows there is a significant difference. What should I do? Am I doing what I am supposed to be doing correctly? In short, I have data collected through many waves (around 5 waves). For your information, the time span between the waves is not consistent. all waves accumulated to 163 respondents (SMEs). What should I do?


title: How to present the Bayesian approach

text: I am writing a dissertation in which I use maximum likelihood estimation and an alternative Bayesian approach. I have written up the maximum likelihood estimation approach. However, I need some advice on presenting the Bayesian part. Can someone suggest me a source which can shape my thoughts so that I can write up the Bayesian part. I need to present the fitted values of a linear model and a logit model and also their Bayesian intervals via plots as the sample size is too large to be given in the dissertation. I want to include winBUGS plots. What plots can be useful to give the fitted values and their Bayesian intervals. Any advice would be greatly appreciated.


title: PCA before random forest classification

text: Does it make sense to do PCA before carrying out a Random Forest Classification? I'm dealing with high dimensional text data, and I want to do feature reduction to help avoid the curse of dimensionality, but don't Random Forests already to some sort of dimension reduction?


title: Creating fuzzy values for binary data

text: I am a little bit rusty at this sort of thing, so I will not be offended if you tell me this is a stupid question. I am doing a logistic regression where my dependent variable is whether or not a person owns a particular product. Among the variables in the model is an indicator of marital status, 1 for married, 0 for not married. This is reliable data taken directly from an application the client filled out, for example. In many cases, the data is missing (maybe the person never filled out an application) but there is a second data source which is somewhat accurate, but not completely so. Would it make sense to fill in the missing values - not with 0 or 1, but something in between? Like if they are probably married, based on the second data source, I give them 0.8 or else a 0.2 if they probably are not married.


title: Split-Plot ANOVA: model comparison tests in R

text: How can I test effects in a Split-Plot ANOVA using suitable model comparisons for use with the `X` and `M` arguments of `anova.mlm()` in R? I'm familiar with `?anova.mlm` and Dalgaard (2007)[1]. Unfortunately it only brushes Split- Plot Designs. Doing this in a fully randomized design with two within-subjects factors:               N  <- 20  # 20 subjects total     P  <- 3   # levels within-factor 1     Q  <- 3   # levels within-factor 2     DV <- matrix(rnorm(N* P*Q), ncol=P*Q)           # random data in wide format     id <- expand.grid(IVw1=gl(P, 1), IVw2=gl(Q, 1)) # intra-subjects layout of data matrix          library(car)        # for Anova()     fitA <- lm(DV ~ 1)  # between-subjects design: here no between factor     resA <- Anova(fitA, idata=id, idesign=~IVw1*IVw2)     summary(resA, multivariate=FALSE, univariate=TRUE)  # all tests ...      The following model comparisons lead to the same results. The restricted model doesn't include the effect in question but all other effects of the same order or lower, the full model adds the effect in question.               anova(fitA, idata=id, M=~IVw1 + IVw2, X=~IVw2, test="Spherical") # IVw1     anova(fitA, idata=id, M=~IVw1 + IVw2, X=~IVw1, test="Spherical") # IVw2     anova(fitA, idata=id, M=~IVw1 + IVw2 + IVw1:IVw2,                           X=~IVw1 + IVw2, test="Spherical")          # IVw1:IVw2      A Split-Splot design with one within and one between-subjects factor:               idB  <- subset(id, IVw2==1, select="IVw1")          # use only first within factor     IVb  <- gl(2, 10, labels=c("A", "B"))               # between-subjects factor     fitB <- lm(DV[ , 1:P] ~ IVb)                        # between-subjects design     resB <- Anova(fitB, idata=idB, idesign=~IVw1)     summary(resB, multivariate=FALSE, univariate=TRUE)  # all tests ...      These are the `anova()` commands to replicate the tests, but I don't know why they work. **Why do the tests of the following model comparisons lead to the same results?**               anova(fitB, idata=idB, X=~1, test="Spherical") # IVw1, IVw1:IVb     anova(fitB, idata=idB, M=~1, test="Spherical") # IVb      Two within-subjects factors and one between-subjects factor:               fitC <- lm(DV ~ IVb)  # between-subjects design     resC <- Anova(fitC, idata=id, idesign=~IVw1*IVw2)     summary(resC, multivariate=FALSE, univariate=TRUE)  # all tests ...      **How do I replicate the results given above with the corresponding model comparisons for use with the`X` and `M` arguments of `anova.mlm()`? What is the logic behind these model comparisons?** EDIT: suncoolsu pointed out that for all practical purposes, data from these designs should be analyzed using mixed models. However, I'd still like to understand how to replicate the results of `summary(Anova())` with `anova.mlm(..., X=?, M=?)`. [1]: Dalgaard, P. 2007. New Functions for Multivariate Analysis. R News, 7(2), 2-7.


title: What is the Bayesian counterpart to a two-sample t-test with unequal variances?

text: I am looking for the bayesian counterpart of the two-sample t-test with unequal variances (the Welch test). I am also looking for a multivariate test, like Hotelling's T statistic. References appreciated. For the multivariate case, suppose that we have $(y_1,\cdots,y_N)$ and $(z_1,\cdots,z_N)$, where $y_i$ (resp $z_i$) is a shortcut for a sample mean, sample standard deviation and number of points. We can assume that the number of points is constant across the whole dataset, the standard deviation the same for all $y_i$ (resp $z_i$) and that the sample means of the $y_i$ (resp $z_i$) are correlated. If you plot the sample means, they follow each other and by connecting them, you get a smooth varying function. Now on some parts the $y$ function agrees with the $z$ function, but on others it doesn't, because $\frac{mean(y_i)-mean(z_i)}{std(y_i)+std(z_i)}$ becomes big. I would like to quantify this statement.


title: Measuring individual player effectiveness in 2-player per team sports

text: I've got a spreadsheet of some team scores. First team to 10 points wins. There are 2 players on each team. The players play with different teammates all the time, although they are not chosen perfectly randomly. No individual scores are kept. So basically we have Bill and Bob beat Andy and Alice 10-4 Jake and Bill beat Joe and John 10-8 ... Is it possible to come up with some ranking for the _individual_ players, based on all the available match data. Basically, to see how much each player contributes to each game in terms of points or relative to the other players? Thanks


title: What article or book clearly states that one can't use protected t-tests in within subjects ANOVA?

text: I have been asked to use protected paired samples t-tests in an analysis. The requester states that if I do not use the overall MSe from my (1 factor with four levels) within-subjects ANOVA when conducting my paired samples t-tests there is not really any protection from the ANOVA. As I remember it, in a between subjects ANOVA this procedure is only defensible if the homogeneity of variance assumption is met. It seems like a likely extension to within-subjects ANOVA might be that this is only permissible if there is no violation of sphericity. Because there are violations in this dataset, I have chosen to apply the Huynh-Feldt correction for sphericity. Regardless, if anything such an approach seems anti- conservative as it provides more degrees of freedom in the denominator. In addition, the help file in R for pairwise.t.test says that "pooling does not generalize to paired tests". The purpose for my planned comparison t-tests is merely to identify the differences between conditions that have resulted in a significant ANOVA. I would like to be able to justify my reasons for rejecting the pooling of error variances, but I am unable to find a citation that clearly states that such an approach is inappropriate. Does anyone know of one? Alternatively, why is my thinking on this issue incorrect?


title: Interpretation of coefficients in polynomial regression for predictive modeling

text: I am building a predictive model (binary target variable) in the financial services industry. One of the (many) potential predictors I am adding to the model is related to the customers checking account balance trend (longitudinal balance). I'd like to capture if the balance is increasing or decreasing and how much. I have access to end of month balances going back a ways. One of the things I was considering is to, for each customer - fit a polynomial regression and include the coefficients into my predictive model. In R, an example of a single customer:               balances <- c(657709,620729,713637,619224,558238,572402,536548,0,0,0)     time <- seq(1:10)     mod <- lm(balances~time+I(time*time))     mod$coefficients[2:3]          mod$coefficients[2:3]     time          I(time * time)      61239.99      -13317.43       **Questions:**   1. Thoughts? Of course the fit can be very poor, but as a global process to include into a predictive model does it have merit? Is there a better way?   2. It seems I have seen description of these coefficients in terms of velocity and acceleration, but I cant find it anywhere. Is this a true interpretation of them?


title: Multiple value Chi Squared test in SPSS

text: I am currently doing a study into factors that affect mortality in a particular type of bacterial infection. I would like to compare certain patient characteristics. One of which is the specialty that the patient is under. I currently have a variable in SPSS that has three values according to the three specialties I have grouped my patients in: 'Haematological', 'surgical' or 'other'. I would like to see if there is any statistical significance between each of these specialties and whether the patient survived or died during their hospital stay. I assume because I am dealing with categorical data I have to do a Chi Square test, however I would like to look at whether each specialty is significant individually not just if there is any significant relationship between 'specialties' and death. Is there any way of getting SPSS to do this, or do I need to break them up individually?


title: How to fight underfitting in deep neural net

text: When I started with ANN I thought I'd have to fight overfitting as the main problem. But in practice I can't even get my NN to pass the 20% error rate barrier. I can't even nearly beat my score on random forest! I'm seeking some very general or not so general advice on what should one do to make his NN start capturing trends in data. For implementing NN I use Theano Stacked Auto Encoder, the code from tutorial that works great (less than 5% error rate) on classifying MNIST dataset. It is a multilayer perceptron, with softmax layer on top with each hidden later being pre-trained as autoencoder (fully described at tutorial, chapter 8). The number of input features are ~50 and output classes ~10. NN has sigmoid neurons and all data normalized to [0,1]. I tried lots of different configurations: number of hidden layers and neurons in them (100->100->100, 60->60->60, 60->30->15, etc,etc), different learning and pre-train rates, etc. And the best thing I can get is %20 error rate on validation set and %40 error rate on test set. On the contrary, when I try to use Random Forest (from scikit-learn) I easily get %12 error rate on validation set and %25(!) on test set. How can that be that my deep NN with pre-training behaves so badly? What should I try?


title: Comparing abundance vs sites

text: I need to apply some basic stats to see if there are significant differences in total abundances of these three communities (X, Y and Z) with duplicate counts at three different locations (A, B and C). Please do suggest best statistical test and how to apply.                         Replicates    Location A                 Location B            Location C     Community X   1              2300                     14380                   1230                   2              2098                     11240                   1370     Community Y   1               830                       300                     86                   2               928                       240                     95     Community Z   1              1870                      1450                    562                   2              2016                      1238                    482      Any help would be very much appreciated!


title: How to visualize the comparison of 3 different types of distances among objects

text: I need to calculate the difference between six time series in three ways:   * Time Series are: Kamel, Dumper, Graben, Traktor, Generator    * Methods are: Euclidean distance, Manhattan distance and maximum distance For example, I show the differences between these time series for Euclidean distance in a table like this: ![enter image description here](http://i.stack.imgur.com/w2uKG.png) (I have two more such tables.) Now I want to visualize these tables in one diagram, with a different color for each method. What data visualization method should I use?


title: Choosing the best set of keywords

text: I have a dataset of tweets collected using twitter streaming API on a particular topic (say 'football') using around 40 keywords. Now if I'm going to track the same topic (football) in future how do I determine the best set of keywords to query. Any researchers working on these areas? How do I solve this keyword ambiguity? What are the sophisticated and state of the art techniques?


title: Interaction in a repeated measures logistic regression

text: I want to test for the effect of interaction between two variables using a repeated measures logistic regression but I don't know how. I am using Stata and the following lines:               xtset participant                 xtgee FIRST_VAR SECOND_VAR, fam(binomial) link(logit)


title: Phylogenetic tree with known historical sequences?

text: I have a multiple sequence alignment that I'm using to construct a phylogenetic tree. Usually, phylogenetic trees are constructed under the assumption that the input sequences are all from the present day -- thus, it forces them all to be leaf nodes. However, my problem is a bit different. I have genome sequences collected at various dates from about 1960 to the present. Because of the extremely rapid rate of evolution (since I'm dealing with viruses), its likely that many of the older sequences are actually ancestors of the more recent ones. Is there any kind of phylogenetic structure that can represent this kind of relationship? Or can I deduce this from a standard phylogenetic tree?


title: How do I express how precise the results of my simulation are?

text: I am working on an UNO Simulator to measure the effectiveness of different strategies on an players win percentage. I would like to know the proper way of expressing the results of my testing. For example, I ran the simulator using two identical random player strategies, with different sample sizes. For each sample size, I ran the simulator 20 times and calculated the Mean and Standard Deviation for each set of 10 thousand - 5 million sample games.                     |  Mean  |Std Dev |           +--------+--------+        10k|0.495390|0.027331|        50k|0.502642|0.013702|       100k|0.500470|0.005940|       200k|0.501116|0.005152|       500k|0.502088|0.003186|     1,000k|0.502720|0.001989|     5,000k|0.502542|0.000790|      I wanted to use these examples to express how confident I am that later tested strategies are precise to some fixed decimal Win:Loss ratio. **Can I simply calculate the CI or margin of error (1.96 * Standard Deviation/SQRT(sample size))?** **Would this mean that I am 95% confident that the true Win:Loss ratio is +/- (sigma/SQRT(n))?** **Is this the proper way of expressing the precision of simulations results, or is their a more accurate way of expressing precision?**


title: Levels or Independent Variables in ANOVA

text: There are four different types of mobile phones (same brand) that are used as stimuli in an experiment.   * In the first example, if "mobile phone" is the independent variable in the study, which one of the two is appropriate? (i) mobile phone is the independent variable with four levels or (ii) this study has 4 independent variables (since there are 4 different mobile phones) The same webpage is modified into two versions: one with a blue background, and the other version with an orange background.   * In the second example, if webpage is the independent variable, can we say that the study has one independent variable with 2 levels or does it have two independent variables? Kindly clarify.


title: Difference between Maximum a posteriori and penalized likelihood

text: Can anyone explain to me the difference between penalised likelihood and maximum a posteriori? I read a paper where the likelihood function is $$L(\theta_1, \theta_2,\theta_3 ; x)=f(x|\theta_1, \theta_2) f(\theta_2|\theta_1,\theta_3)f(\theta_3)$$ or $$\ell(\theta_1, \theta_2,\theta_3 ; x)=\log(f(x|\theta_1, \theta_2)) +\log(f(\theta_2|\theta_1,\theta_3))+\log(f(\theta_3))$$ and the authors say they use maximum penalised likelihood to find $\widehat\Theta$. I thought that would be maximum a posteriori and penalised likelihood would use a likelihood function of the form. $$\ell(\Theta ; x)=\log(f(x|\Theta)) - \lambda g(\Theta)$$ Any thoughts?


title: How many users will connect to my Server?

text: I've implemented a backend which expose a rest API and it's now running on my own server. I know it will be used by 5 groups of 20 users. Each group will use the server for 2 hours a week. If a group starts working, all the users will work at the same time, and for the following 2 hours, on the server. The choice of working hours is indipendent between groups. Each week is composed by 5 days of 6 hours per day. Can I estimate the probability of having N users connected at the same time, **during the working hours**? Can you help me formalizing this problem? The objective of this model is to estimate which is the maximum number of users should I expect, at which probability. I intuitively understand that having in a given instant 100 users connected is less frequent than having 20. But "how much"?


title: .632+ bootstrap estimator-defining gamma in continuous variable case

text: I am trying to implement the .632+ bootstrap estimator for internal validation as proposed by Efron and Tibshiraini 1997. Looking at the paper, I can see how gamma is defined in the case of classification problem (equations 26 and 27), but is there an equivalent that is typically applied when predicting a continuous variable? Is the .632+ estimator applied in the context of predicting continuous variables? Thanks


title: Is sample homogeneity an assumption of regression analysis?

text: I have assumed (i.e. I think I was taught, longer ago than I can remember) that regression analyses assume that a sample is homogeneous. If it is not, then the appropriate thing to do is either add dummy variables to code for the different groups included in the sample, or carry out an ANCOVA to test whether group parameters are equal. Does ignoring the heterogeneity of a sample invalidate a regression analysis?


title: What's the difference between the output of Ridge regression with linearRidge in ridge package and lm.ridge in MASS package?

text: When I fit a ridge regression, in the example that I show in below, with the same lambda=0.1 the estimated coefficients from the ridge package and MASS package are different. What's the difference? for example when I want to fit dependent variable z to independent x variables as I show below:               x<-rnorm(50,mean=5,sd=1)     z<-matrix(ncol=10,nrow=50)     for (i in 1:50)     {     z[i,]<-rbinom(10,1,.7)     }     massridgefit<-lm.ridge(x~z,lambda=.1)     ridgefit<-linearRidge(x~z,lambda=.1)      when we estimate coefficients with:               coef(ridgefit)     coef(massridgefit)       the results is:               (Intercept)          z1          z2          z3          z4          z5      5.23671868  0.08236113  0.14447491 -0.27670037 -0.66643958 -0.02028575           z6          z7          z8          z9         z10      0.01164812  0.05839092 -0.46193876  0.26368953  0.48542056       and                               z1          z2          z3          z4          z5      5.26038037  0.11131652  0.18158273 -0.30602658 -0.76933007 -0.03197175           z6          z7          z8          z9         z10      0.01836979  0.07302220 -0.54970233  0.29167218  0.56262584      respectively. What's the difference?


title: Is this Markov chain irreducible?

text: I am trying to figure out if this Markov chain is irreducible and if it is aperiodic and why or why not. ![enter image description here](http://i.stack.imgur.com/i70tF.png) For me it is not irreducible Markov chain because you cannot go from state 2, 5 or 6 to any of the states 1, 3, 4. I don't know if my argument is good.


title: Ordinal logistic regression with likert scales

text: I'm currently have a bit of difficulty determining how to analyze this data via logistic regression analysis.                - Q18 = DV (satisfaction score ranging from 1-10)      - Q10_1 = IV (Customer Service likert score from 1-5)      - Q10_2 = IV (Sales likert score from 1-5)      - Q10_3 = IV (Performance likert score from 1-5)      - Q10_4 = IV (price likert score from 1-5)       - Q10_5 = IV (proposal likret score from 1-5)      - Q10_6 = IV (collateral likert score from 1-5)      - Q10_7 = IV (reporting likert score from 1-5)      - Q10_8 = IV (manager likert score from 1-5)      My guess is that you need to use an ordered logistic regression model but i'm not sure what to factor in my formula just the DV or everything? Which equation is correct here? > nps.olr <- polr(data = cs_aggmean, formula = factor(Q18) ~ Q10_1 + Q10_2 + > Q10_3 + Q10_4 + Q10_5 + Q10_6 + Q10_7 + Q10_8) or > nps.olr <- polr(data = cs_aggmean, formula = factor(Q18) ~ factor(Q10_1) + > factor(Q10_2) + factor(Q10_3) + factor(Q10_4) + factor(Q10_5) + > factor(Q10_6) + factor(Q10_7) + factor(Q10_8)) After that I'm have trouble interpreting the results. For the first model the odds ratio I believe is this after I exponentiate the coefficients: > exp(nps.olr$coefficients)   >  Q10_1 = 1.834354   >  Q10_2 = 1.354964   >  Q10_3 = 3.259454   >  Q10_4 = 1.269431   >  Q10_5 = 1.326062   >  Q10_6 = 1.432196   >  Q10_7 = 1.424732   >  Q10_8 = 1.010827   > I appreciate any guidance here and of course just let me know if I need to supply more information! I should mention that I'm using R for software and that I'm less interested in making a predictive model and more in making recommendations on how to increase satisfaction from these variables.


title: Standard deviation of normalized data

text: I have a data set $y_i$ (where the $y_i$ are photon counts in time period $i$, $i=1,2,...,n$, assumed Poisson), with an estimated standard error $s_i$ (= $\sqrt{y_i}$) for each count. For some reason, I normalize the data set by the average count rate $\bar y$. $\bar y = \frac{\sum y_i}{N}$ where $y_i$ has $i=1..N$ The standard deviation of $\bar y$ is calculated as $\sqrt{\frac{\sum (y_i - \bar y)^2}{N-1}}$ To normalize I simply do: $z_i = \frac{y_i}{\bar y}$ At this point how do I calculate the standard error for the $z_i$? Is error propagation the right thing? $\text{se}(z_i) = z_i \times \sqrt{(\frac{s_i}{y_i})^2 + (\frac{\text{se}(\bar y)}{\bar y})^2}$ It seems to me that the error stay of the same order as before, while the count rates increase their value about one magnitude. So I am probably missing something. So, if I print the ratio of the count rate values divided by their standard deviation, before the normalization, I obtain:               [  1.06904497   1.07222193  -0.306786     2.22599555   5.04049535     10.44367859  10.37041246   4.71728177  10.85418506  10.99314159     10.18889392   9.20449287  13.24244513  10.70825227  15.64406957     15.41271307  10.06729494  12.88619079  14.76416192   3.90199486     13.33680974  11.57703491   9.38122633  11.53686373  13.32397254     -0.78308901]      And after the normalization (plus error propagation):               [  2.27976751e+00   2.29330204e+00   1.88194867e-01   9.79928633e+00     4.80280889e+01   1.74659904e+02   1.72699185e+02   4.23542991e+01     1.85693431e+02   1.89443907e+02   1.67854381e+02   1.41988475e+02     2.50472302e+02   1.81762548e+02   3.14047817e+02   3.08069283e+02     1.64619984e+02   2.40831046e+02   2.91120670e+02   2.94284485e+01     2.53020286e+02   2.05263389e+02   1.46572585e+02   2.04172711e+02     2.52673822e+02   1.69702062e+02   2.48534464e+01   1.65477547e+01     4.12770323e+00   4.47260177e-01   1.22474257e+00]      Any suggestion really appreciated. PS: Sorry if the formulas are not well displayed, I don't know how to set the "formula" environment.


title: Problem when creating matrix of values based on covariance matrix

text: I want to simulate a data set with similar covariance structure as my observed data (which is a SNP by gene p-value matrix, dim ~600k*8368), and have calculated a covariance matrix (dimensions 8368*8368). So far I've tried two approaches to simulating data (all in R): `rmvnorm` from the mvtnorm package, and by using the Cholesky decomposition. The problem is that the resulting covariance structure in my simulated data is very different from the original supplied covariance vector (but the code used is accurate as far as I can tell). Lets just look at some of the values:                   > cov8[1:4,1:4] # covariance of simulated data             X1          X2         X3         X4         X1 34515296.00    99956.69   369538.1  1749086.6         X2    99956.69 34515296.00  2145289.9  -624961.1         X3   369538.08  2145289.93 34515296.0  -163716.5         X4  1749086.62  -624961.09  -163716.5 34515296.0         > CEUcovar[1:4,1:4]              [,1]         [,2]          [,3]         [,4]         [1,] 0.1873402987  0.001837229  0.0009009272  0.010324521         [2,] 0.0018372286  0.188665853  0.0124216535 -0.001755035         [3,] 0.0009009272  0.012421654  0.1867835412 -0.000142395         [4,] 0.0103245214 -0.001755035 -0.0001423950  0.192883488      So the distribution of the observed covariance is very narrow compared to the simulated data. None of the eigenvalues of the observed covariance matrix are negative, and it appears to be a positive definite matrix. Here is what I did to create the simulated data (very straightforward):                   Chol <- chol(CEUcovar)         Z <- matrix(rnorm(20351 * 8368), 8368)         X <- t(Chol) %*% Z         sample8 <- data.frame(as.matrix(t(X)))         > dim(sample8)         [1] 20351  8368         cov8=cov(sample8,method='spearman')      (earlier I've also tried: `sample8 <\- rmvnorm(1000,mean=rep(0,ncol(CEUcovar)), sigma=CEUcovar, method="eigen")` with as 'bad' results, much larger covariance values in the simulated data)   * Any ideas of WHY the simulated data have such a different covariance?    * Any experience with similar issues? Would be happy to supply the covariance matrix if anyone wants to give it a try.    * Is this an impossible task statistically?   * Any suggestions? Anything apparent that I left our or neglected? I am aware that my supplied covariance matrix might differ from the true underlying covariance of the population. But my goal is to simulate a random data that captures (as much as possible) of the covariance that this observed data has, no matter how weird that covariance is. _Clarification - my supplied covariance vector has values in [-10^-2, 0.3] while the resulting data has covariance of [-800000,5000000]. This is not what I am expecting. I wonder why this is the case, and how this might be explained._ Any advice would be highly appreciated.


title: Statistical significance between two groups

text: Heart                  Group-1                         Group-2                    Case1    Case2   Case3          Case1    Case2   Case3     İnflammation    -       -       -              +++      ++      ++     Hemorrhage      +       -       -              ++       +++     ++     Fibrosis        +       +       -              ++       +++     ++     Necrosis        ++      ++      +              ++       +++     ++          (-;absent, +;mild, ++;moderate, +++;severe)      A sample table is seen above. I have 3 cases in each groups. How can I calculate statistical significance between two groups. Should I calculate statistical significance for each lesion?


title: Interpreting t-value as a measure of evidence

text: 1. Came across an article on the web on "evaluation of training effectiveness". The author suggests that the "t-value" obtained from a "Paired t-test" conducted using pre-test and post test scores can be used to quantify the effectiveness of training.    2. The t-value was termed as the "Index of Learning" and the author claimed the following:   a. An "Index of learning" (t-value) between 0 to 1.5 indicates no evidence of learning.   b. A t-value between 1.5 to 2.0 indicates some evidence of learning.   c. A t-value between 2 to 3 indicates strong evidence of learning.   d. A t-value above 3 indicates very strong strong evidence of learning.   3. In my limited understanding of the subject, the t-value, critical t-value and p-value of a t-test would depend on the significance level and degrees of freedom of the associated data. I have thus not understood the rationale behind the generalization made by the author (Para 2 (a) to 2 (c) above). I would be extremely grateful for any explanations. Thanks.


title: Can machine learning find all sort of crazy connections?

text: If you try a real thoroughly won't a computer find all sort of silly patterns?   * Messages from ETs in the bible   * rainy Sundays in China or Australia -> the chances of your sport team win   * reading many novels -> having a gay son If you have a huge (and growing amount) of information, would things start to match in a meaningless way? What safe-guards do you put in place to get rid of these associations?


title: Why is volatility an important topic in financial econometrics?

text: I do not know if it is totally off-topic, but I thought it might be useful to have opinions and an aggregate answer about why volatility is an important topic in financial econometrics. I think it started with portfolio theory and the need to understand the properties of the underlying second moment of the asset returns. Subsequently the Black-Scholes formula and the popularity of derivatives made this entity very important in Finance.


title: semisupervised classification training on all or part of unlabeled data

text: I have 3 sets of data.   1. A positively labeled dataset.    2. An unlabeled dataset that has for sure positive (around 75%) and negative data.   3. An unlabeled dataset that has for sure positive data and maybe negative data. I'm interested in finding negative observations in this last dataset. There will be no future new (test) data. Im trying out some semisupervised logistic regression techniques using shift parameters (http://arxiv.org/abs/1305.4987).   1. Since I wont have future data should I include all my samples in building my classifier and classify my training data to determine the labels of dataset 2? Or should I train my classifier on only dataset 1 and 2 and then determine the labels of dataset 3?   2. If I include all data would I still have a problem with overfitting? Do I still need L1 regularisation or something?    3. Is my choise of training algorithm appropiate for this problem?   4. Im also strugeling what a good performance measure would be. Missclassification error is probably not ideal. best regards


title: How to specify a non MLE prior for a naiveBayes model in R?

text: I am just learning R, and am aware that the package e1071 has a `naiveBayes` method that takes in predictor and class membership, and estimates the class prior using the relative frequency (ML estimate). I would like to see the impact of having a sample that is not balanced across classes (let's say male, female are represented in 3:7 ratio in the available data), when I know that for the population I am interested in, the prior of male, female would be much closer to 0.5:0.5 than 0.3:0.7. The skewed representation is due to practical problems in sampling/surveying. I've tried doing a down-sampling of the female subset to make the training set reflect the 1:1 ratio. Is there a way of using `naiveBayes` with a pre- specified prior that is not reflected in the relative frequency in the data? I've looked through the documentation on `e1071::naiveBayes` and there does not seem to be any way of specifying my own prior. Is there another package that does that?


title: weighted conditional expectation in adaboost

text: I'm looking at this paper, http://www.stanford.edu/~hastie/Papers/AdditiveLogisticRegression/alr.pdf around pages 10-11 (marked 346,347 in the pdf). This notation is introduced $$ E_w[g(x,y)|x] := \frac{E[w(x,y)g(x,y)|x]}{E[w(x,y)|x]} $$ where originally, we are attempting to minimize this: $$ E[w(x,y)g(x,y;f,c)] $$ with respect to f and c where $w(x,y)=\exp(-yF(x))$, $g(x,y)=1-cf(x)+c^2/2$. Earlier, it was shown in minimizing $E[w(x,y; F(x))]$ with respect to F that it is sufficient to minimize the conditional, $E[w(x,y;F(x))|x]$. Fine. This makes sense so far, but now in minimizing $E[w(x,y)g(x,y;f,c)]$ it seems as though we could take the same route and minimize $E[w(x,y)g(x,y;f,c) | x]$. This would be equivalent to $E_w$ without the denominator. I see that scaling should not affect the minimum found since it does not involve $f$ or $c$. I also see that this scaling gives $E_w$ the usual linear expectation properties. Is this introduced only for convenience or is there some other justification?


title: PCA on standardized variables: how do eigenvectors relate to original, unstandardized variables?

text: PCA performs a linear transformation on a data set to obtain a new data set, this time with eigenvector basis and eigenvalue loadings. If Z is our data set, P our linear transformation, and Y our new data set with eigenvector basis, the transformation is: P*Z = Y . but, quiet often Z is not our original data set. If the units of the variables are different, we must perform some sort of standardization/normalization transformation, say K, so that they are comparable. If K is our transformation matrix, then to obtain the new, standardized data set, we perform: K*X = Z My question is, what are the loadings in terms of our original data set X? In more practical terms, if b1 is the first measurement type of Z, PC1 might look like.. b1 5.0   b2 2.7   b3 3.0   ...   bn 10.0 Where the values of b1->bn are given by P. But what is PC1 in terms of the original data set X, with measurements a1->an?    a1 ?   a2 ?   a3 ?   ...   an ? These should come from some new linear transformation P_new. I would assume we can use basic matrix algebra to figure this out, but I want to make sure I'm 100% correct here: if K*X = Z, and P_old*Z = Y, we're searching for: P_new*X = Y we should be able to use... P_new*X = Y or P_new = Y*X^-1 ? Is this correct? What if X is not square (and it most certainly is not, there are almost always more trials than variables in any good data set)? Will computational solutions make sense? I've searched through 23/34 pages on PCA, and was not able to find an answer to this question. If its out there, my apologies, I tried.


title: Updating ARIMA models at frequent intervals

text: I am using an ARIMA model to create a model for correlated errors from my regression model. I am using the `auto.arima` function from the forecast package in R. I am able to get more data at some frequent interval after the regression model is created, so I get more values for the correlated errors. My question is how do I update the ARIMA model with a gap in time interval between readings.


title: A probability question that needs some clarification

text: I have a probability question (if it is too basic, my apologies since I am learning), I can not get my head around. Suppose I concatenate the sentences from many books and shuffle the sentences (so we do not know the book that a sentence is from) to have a "sentence database". We assume there is no duplicated sentences, and we knew the books are written by two different writers. Now, given a phrase (sequence of words), we can find this phrase at 4 sentences in our database. And we are also given that for these 4 sentences:               sentence 1 has 20% probability is written by writer 1          sentence 2 has 60% probability is written by writer 1          sentence 3 has 50% probability is written by writer 2          sentence 4 has 40% probability is written by writer 2      The question is what is the probability of the given phrase is written by writer 1 ~~and~~ **or** writer 2? My calculation is for the probability that the phrase is written by writer 1: $$\frac{20\% + 60\%}{20\% + 60\% + 50\% + 40\%} = 0.470588235$$ and for writer 2: $$1 - 0.470588235 = 0.529411765$$ My intuition tells me I am doing it wrong. But i do not where is wrong, could you please point it out. Edit: I thought this question is similar with calculating the probability of a outcome when flipping a coin. So if I flip a coin 4 times, with outputs H T T H, then the probability of H is 1/2. So the formula is $$\frac{event\ occurred}{total\ number\ of\ trials}$$ Edit 2: _Where do the probabilities of the sentences come from?_ It is coming from a machine learning algorithm (SVM). _What do they mean?_ Because the sentences are shuffled, we lost the information about which writer wrote this sentence. However, we know the (sentence) writing styles of these two writers, so we trained a SVM by using the known writing styles as features to give each sentence a probability, the probability means this sentence has x% probability is written by writer y. As the task is getting more difficult now. We are just given a phrase to judge (I mean assign a probability to this phrase, in case you think judge might means binary decision) which writer wrote it. (You might say can you not train another SVM to get the probability, just like what you did for sentence probabilities. Assume we can not do it.) I had modified my question, it should be "The question is what is the probability of the given phrase is written by writer 1 **or** writer 2?"


title: Delete observations for regression modelling

text: I am implementing robust regression on some data and when i tested the accuracy on a hold set of data, the accuracy was very bad. I have about 500 000 observations for the regression model and i deleted about 20% of the observations in such a way that each observation should include at least 3 of the 10 independent variables. What i realised, if i delete the observations this way i can improve my accuracy alot. The significance of the variables also stays the same whether i delete the observations or not and the hold out data set contains observations that only has one or two of the independent variables. But i would like to know whether there is a formal method for deleting observations?


title: Metropolis Hastings when acceptance rate is not a probability

text: I need to implement a Metropolis Hastings where the acceptance probability $\alpha$ is not a probability but a logarithm of a score. The logarithm of a score is a negative float. In original algorithm i have:               alpha = compute acceptance probability with proposal distribution     Compute r = min(1,alpha)     Sample u in U(0,1)     Set new sample to         x' if u < r         x if u >= r      min() between 1 and alpha works fine with probability but if alpha is like -30 -(-54) or -23 - (-10), in the first case it's ok (i accept the new state), but in second case ? Observation: I use the logarithm of a score then my acceptance probability is $$\frac{\text{score(new_state)}}{\text{score(old_state)}}$$ and with logarithm is $$\ln\text{(score(new_state))} - \ln\text{(score(old_state))}$$


title: analysis of alpha lattice design

text: Are there R some packages for analysis of alpha-lattice (0,1) designed experiment? My experiment includes two factors with two levels. The tests are performed in two replications with twenty blocks. Every block consists of 10 plots.


title: If betas of two stocks are exactly the same

text: This is a very basic question - If two stocks have the same beta over same time period, does it mean they are 100% correlated over that time period? Beta, in a CAPM framework is defined as Cov (R1, M) * Var (M). Where R1 is the return vector of a security 1 and M is the market return vector. Equating two betas mean Corr(M, R1) * Std (R1) = Corr (M, R2) * Std (R2). I'm not really sure where to go from here - the standard deviations of R1 and R2 might not be equal, and I'm not sure what the relation, if any is between the Corr (M, R2) and Corr (M, R1). According to the paper below, correlation is not transitive. If R1 and M and perfectly correlated, and R2 and M are perfectly correlated, it doesn't necessarily mean R1 and R2 are perfectly correlated. Any thoughts? http://www.jstor.org/discover/10.2307/2685695?uid=3739832&uid=2&uid=4&uid=3739256&sid=21102787273757


title: How to predict binary outcome from a glmm model

text: Suppose I fit a generalized mixed logistic model such like that:               set.seed(2014)     require(lme4)     df<-data.frame(id=rep(1:5, c(8,10,12,14,15)),                     out=c(rbinom(8,1,0.1), rbinom(10,1,0.3),rbinom(12,1,0.1),rbinom(14,1,0.05),rbinom(15,1,0.1)),                    age=rnorm(59,50,10),                    gender=rbinom(59,1,0.5))     fit<-glmer(out~age+gender+(1|id),data=df,binomial)     df$predicted<-predict(fit,type="response")     df$pred.binary<-with(df,ifelse(predicted>=0.5,1,0))     apply(df[,c(2,6)],2,sum)             out pred.binary                9           0       What's the best way to predict a binary outcomes by accounting the random effects and with an optimized threshold based on the original data set?


title: Adding 2nd level variable into Multi-level Modelling in Stata

text: I'm used to HLM 7 software and now I'd like to switch to Stata, for multilevel modelling (xtmixed). To give an example imagine I have students(level1) nested within schools (level2). In HLM I can easily add a second level variable (for instance schools beauty) choosing the equation of a coefficient (or intercept) which is at level 2. Imagine I want to include the effect of school beauty in B0 (intercept). How can I do that in Stata?


title: Predicting time to finish

text: Out of curiosity, I want to understand how to model this problem. I've been hearing people suggest the use of linear regression but **I am not sure how to encode this problem** (included my attempt below) in R as I am a complete beginner in this area. I have a task that can be done any number of times (each individual instance is a task instance). Everytime the task completes 1%, I recorded the time elapsed since the task's start time. Therefore, for each task, I will have 100 points (100 1% increments) at which I recorded the time elapsed. Given that I have this data for many instances, is it possible to predict the finish time for this task when a new task instance is given?                     TaskID Percent TimeElapsed        1:      1       0   0.2035333        2:      1       1   0.2062833        3:      1       2   0.2137167        4:      1       3   0.2180833        5:      1       4   0.2490833       ---                                3127:     31      96   4.9391667     3128:     31      97   4.9970500     3129:     31      98   5.5644500     3130:     31      99   5.6532667     3131:     31     100   5.8359833      A quick look at the task behavior (below) tells me there is a bit of a variance in how the task behaves so its hinting that the output should not just be a time prediction but rather a time prediction with some confidence? In addition, I'm thinking just using the information about the current progress of the task might not be sufficient - the task may have slowed down in some its previous progress points so the finsh time would be affected. Therefore, this information should somehow be encoded into the model? ![enter image description here](http://i.stack.imgur.com/eiEKh.png) I am particularly interested in how to do this using R. I included my initial attempt at using linear regression here but the result does not look good to me. Any suggestions on how to improve this or use some other methods? I have given the output of dput (on a data table: `install.packages("data.table")`) on pastebin. If you want a data.frame instead, please see this paste instead. **EDIT:** Attempt at using linear regression The thick black line is the median at every point. The thick red line is the regression line fit to the median line. ![enter image description here](http://i.stack.imgur.com/mZFZd.png)


title: Paired t-test and random intercept model

text: It's known that the paired t-test is equivalent to the random intercept model in the sense that the data for paired t-test is actually two repeated measures on the same subject and we can use a indicator variable as a covariate in the r.e. model stands for the two measures and with the random intercept to explain the within subject correlation (repeated measures). Can anyone give a formal proof of this and provide a simple data example? Thanks.


title: Big picture on survival analysis and life data analysis

text: I have heard of survival analysis and life data analysis, but don't quite get the big picture. I was wondering what topics they are covering? Is it pure statistics, or just application of statistics on some specific area? Is life date analysis part of survival analysis? Thanks and regards!


title: Rule of Thumb for minimum length of time series for Autocorrection estimation

text: I had a related question answered here: Rule of Thumb for minimum length of time series for AR(1) estimation However the answer gives rise to a new question. I want to be able to estimate the Auto correlation at lag 1 for different step sizes. While I agree that in theory if the lag 1 auto correlation for step size one is $\phi$ then the lag 1 auto correlation for a step size h should be $\phi^h$. However I want to test this empirically to see if it is in fact the case. From the testing I have done it seems that for smaller step sizes the theoretical relationship holds in the empirical data set, but as the step size gets larger it starts to deviate (and also to oscillate) so I think at least partly what I am seeing is the noise associated with estimating the lag 1 auto correlation on a small time series. Is there a formula or even a rule of thumb that I can use to decide what is an acceptable time series length for a lag 1 auto correlation estimation? Thanks Baz


title: k-medoids algorithm with incomplete distance matrix

text: I want to apply k-medoids algorithm using an incomplete distance matrix as input. How can I handle the lack of information of this matrix? Just ignoring the missing distances? Or is there a better way?


title: Statistic (Linear normal model): $X_{hi} \sim N(\alpha_h+ \beta t_{hi}, \sigma^2)$. How to calculate $C_{95}(\beta)$ and $t$-test for $\beta = 2.5$?

text: > Statistic (Linear normal model): $X_{hi} \sim N(\alpha_h+ \beta t_{hi}, > \sigma^2)$. How to calculate $C_{95}(\beta)$ and $t$-test for $\beta = 2.5$? I am in the statistical model $X_{hi} \sim N(\alpha_h + \beta t_{hi}, \sigma^2)$ $h=1,..,3, i=1,..,n_h$ (Linear normal model) where my observations come from three models of linear regression. I have already tested for equal variance and beta value. I'm asked to calculate a $95$%-confidence interval for $\beta$ and to test $\beta = 2.5$. How do I compute these? If I knew the distribution it would be a start.


title: How to calculate Mid-p confidence intervals for SMR, SIR (observed/expected)

text: I need help calculating the Mid-P of the confidence limits for the standardized mortality ratio, or any observed/expected ratio. I am familiar with the Mid-P for a binomial calculation but I am trying to understand how to calculate confidence intervals for a standardized mortality ratio under a Poisson distribution. I've found calculators with formulas such as OpenEpi but I haven't been able to replicate their results. Does anyone know of a source describing this calculation, or of a simple example that could help?


title: Calculating standard deviation after spatially and temporally averaging data

text: I have a spatio-temporal data set with (n,m) spatial and k temporal dimensions. My initial analysis consisted of spatially averaging the data and looking at the time dependent behavior. This resulted in a vector of length k for the mean and standard deviations. Now I want to average the those means, to get an overall scalar value. How do I calculated the standard deviation of this new quantity? Should I be using the pooled variance to calculate it? Or is there another approach that is better suited? Here is an example of the data after spatially averaging. The light regions represent 95% confidence intervals. What I would like to do, is average the last 10 seconds of data (50--60 s) and calculate the corresponding standard deviation so that I can get confidence estimates for the mean. ![](http://i.stack.imgur.com/JvMVa.png)


title: Negative Adjusted R Square is a "good" thing?

text: Working on dissertation. Theory: in K-12 education putting more administrative authority in the state board of education is "better" that leaving it to the local boards. DV: data from 43 states on % kids graduating high school within 4 years. IVs: I have 37 different measures for the types of administrative authority- who is responsible for buying school books, setting school stop/start dates, etc.: 3 (state has complete control), 2 (shared/split) and 1 (locality has complete control). So I fire up SPSS, data from 43 states on % kids graduating high school within 4 years as my DV, plunk in the 37 IVs, use "Enter" as my method (I've been told stepwise is evil, evil, evil) and... > Model Summary > > R = .853 > > R Square = .727 > > Adjusted R Square = -.041 > > Std. Error of the Estimate = .148600607125323 This might be "good" if it means the predictors are useless and the hypothesis is wrong (or at least cannot be proven with this data). It is "bad" if I am getting this because my model stinks. How can I determine which?


title: Estimating entropy of multidimensional variable through dimension reduction

text: Is there anything inherently wrong with trying to estimate the entropy of a multidimensional random variable by first transforming it (by some method) into a single-dimensional variable?


title: Kruskal-Wallis?

text: I have three sets of data (1 set before treatment, second during treatment and the third after), each of which consists of hundreds of trials on a specific task that involves choosing 1 from 7 possible choices, the answer for which may be correct or not (binary categorization). I'd like to determine whether the proportion of correct choices is significantly different as a result of treatment. Which statistical test could potentially help?


title: SAS Scoring on new data

text: I have a training dataset with (x1,x2,x3,y) and these contain some missing observations. I ran proc logistic on this training data and got the parameter estimates (bo, b1, b2, b3). I tried to apply these estimates(bo, b1, b2, b3) on the test dataset which does not contain any missing observation using **proc logistic inmodel** the results show lot of missing data and this is the message from the log _" Work.Test data set are not scored because they have class levels that are missing or are not present in the analysis data set"_ Any suggestion how to score on test dataset without producing missing data even when training dataset contains missing data ? Thanks in advance.


title: Accuracy rate in naive Bayes classification

text: I am trying to use a naive Bayes classification technique to predict fraudsters (`Caller`). My training set of 138 instances has 5 columns viz. `Morning`, `Afternoon`, `Evening`, `Night` and `Caller`. `Morning` has 8 names; the rest all have 3. The names are unique to each column.               > levels(mlcallers$Morning)         [1] "Kelly"   "Larry"   "Mark"    "Nancy"   "Olga"    "Peter"   "Quentin" "Robert"          > levels(mlcallers$Afternoon)     [1] "George" "Harry"  "John"       > levels(mlcallers$Evening)         [1] "David" "Emily" "Frank"         > levels(mlcallers$Night)     [1] "Alex"  "Beth"  "Clark"     > levels(mlcallers$Caller)     [1] "Sally"    "Vince"    "Virginia"          ## Training data set     > summary(mlcallers)     Morning    Afternoon   Evening     Night         Caller        Olga   :29   George:31   David:35   Alex :43   Sally   :47        Peter  :29   Harry :49   Emily:44   Beth :52   Vince   :43        Quentin:22   John  :58   Frank:59   Clark:43   Virginia:48        Robert :13                                                        Mark   :12                                                        Nancy  :12                                                        (Other):21            ## Test data set     > summary(testdata)      Morning   Afternoon  Evening    Night        Kelly  :2   George:7   David:7   Alex :6        Larry  :1   Harry :3   Emily:1   Beth :2        Mark   :1   John  :5   Frank:7   Clark:7        Nancy  :1                                       Olga   :1                                       Quentin:4                                       Robert :5       I need to predict the `Caller` (probable fraudster) in the test data set and also report their confidence. My attempt is as follows:               > library(e1071)     > model <- naiveBayes(Caller~., data=mlcallers)     > predict(model, testdata)     [1] Sally    Sally    Sally    Vince    Sally    Vince    Vince    Virginia    Virginia     Virginia    Virginia    Sally        [13] Sally    Virginia Sally        Levels: Sally Vince Virginia     > predict(model, testdata, type="raw")            Sally       Vince    Virginia      [1,] 0.81806260 0.155576135 0.026361264      [2,] 0.93405930 0.057871797 0.008068906      [3,] 0.82235738 0.122064542 0.055578078      [4,] 0.40028059 0.540385630 0.059333784      [5,] 0.74235622 0.068005993 0.189637784      [6,] 0.14954988 0.723948800 0.126501321      [7,] 0.12730200 0.657333452 0.215364548      [8,] 0.12960601 0.001299552 0.869094437      [9,] 0.02378972 0.001962829 0.974247449     [10,] 0.01420016 0.171655799 0.814144039     [11,] 0.01588595 0.101070603 0.883043443     [12,] 0.82235738 0.122064542 0.055578078     [13,] 0.93569932 0.052176068 0.012124610     [14,] 0.02152574 0.361402108 0.617072156     [15,] 0.74235622 0.068005993 0.189637784      Since the test data has only 15 instances, to find out the accuracy, I have used:               > pred <- predict(model, newdata=testdata, laplace=3)     > table(pred, mlcallers[65:79,5])          pred       Sally Vince Virginia       Sally        1     4        2       Vince        0     3        0       Virginia     0     0        5      ....which gives the highest accuracy in some random trials. My question has 2 parts:   1. Is this the correct approach to the problem?   2. Is there any way I can find out the highest accuracy without me having to randomly select 15 rows from the training data (different random selection like [110:124], yields different accuracy results?


title: How to present a one-way ANOVA when one of the group's variance is "unequal."

text: I'm writing a research paper in entomology that compares 6 categorical variables (species) in terms of one dependent anatomic continuous variable. My samples look very normal by box-plots, etc. Looking at the "largest standard deviation can't be more than twice the smallest standard deviation rule of thumb, one of the six samples has a standard deviation that is less than half the largest $0.01847$ vs. $0.05181$). The other st. deviations are all "OK". Sample sizes range from 13 for several species to 32 for another. Doing a one-way ANOVA on the six samples and ignoring the one species too small st. dev. I get $F > 40$ and a $p$ value of $9.6 \times 10^{-24}$, so I am sure the result is significant. Can I just do the ANOVA anyway. What comment should I make about the one "smallish" variance? Is it a problem to have very unequal sample sizes? (Mine are 25, 22, 13, 13, 13, 32.)


title: Parameters with categorical and gamma distributions in posterior distribution

text: I'm following a very good IPython notebook (the whole list can be found here) in which some sampling techniques are explained. However, I don't understand the use of a categorical variable in a change point model. The posterior distribution for the model is the following: \begin{aligned} P( \lambda_1, \lambda_2, \tau | \mathbf{y} ) &\propto \left[\prod_{t=1851}^{\tau} \text{Poi}(y_t|\lambda_1) \prod_{t=\tau+1}^{1962} \text{Poi}(y_t|\lambda_2) \right] \text{Gamma}(\lambda_1|\alpha,\beta) \text{Gamma}(\lambda_2|\alpha, \beta) \\ &\propto \left[\prod_{t=1851}^{\tau} e^{-\lambda_1}\lambda_1^{y_t} \prod_{t=\tau+1}^{1962} e^{-\lambda_2} \lambda_2^{y_t} \right] \lambda_1^{\alpha-1} e^{-\beta\lambda_1} \lambda_2^{\alpha-1} e^{-\beta\lambda_2} \\ &\propto \lambda_1^{\sum_{t=1851}^{\tau} y_t +\alpha-1} e^{-(\beta+\tau)\lambda_1} \lambda_2^{\sum_{t=\tau+1}^{1962} y_i + \alpha-1} e^{-\beta\lambda_2} \end{aligned} After this, the parameters are given by: $$\lambda_1 \sim \text{Gamma}(\sum_{t=1851}^{\tau} y_t +\alpha, \tau+\beta)$$ $$\lambda_2 \sim \text{Gamma}(\sum_{t=\tau+1}^{1962} y_i + \alpha, 1962-\tau+\beta)$$ $$\tau \sim \text{Categorical}\left( \frac{\lambda_1^{\sum_{t=1851}^{\tau} y_t +\alpha-1} e^{-(\beta+\tau)\lambda_1} \lambda_2^{\sum_{t=\tau+1}^{1962} y_i + \alpha-1} e^{-\beta\lambda_2}}{\sum_{k=1851}^{1962} \lambda_1^{\sum_{t=1851}^{\tau} y_t +\alpha-1} e^{-(\beta+\tau)\lambda_1} \lambda_2^{\sum_{t=\tau+1}^{1962} y_i + \alpha-1} e^{-\beta\lambda_2}} \right)$$ How can we obtain the distributions for these parameters? By marginalizing the posterior distribution? In the case of $\lambda_{1}$ and $\lambda_{2}$ each parameter is described by a Gamma distribution multiplied by some constant, so maybe that is the idea. In the case of $\tau$ I can't see why we obtain a categorical distribution with that argument. Thanks a lot


title: Testing for equality of differences in proportions

text: Say I have the following college enrollment data: **California** In 1975, 188,234 men and 156,612 women were enrolled in California colleges.   In 2005, 194,416 men and 201,334 women were enrolled in California colleges. **Texas** In 1975, 132,261 men and 98,554 women were enrolled in Texas colleges.   In 2005, 140,082 men and 153,305 women were enrolled in Texas colleges. Consider, for each state, the proportion of those enrolled in colleges who were female. How would I go about computing the p-value for the null hypothesis that the 1975–2005 increase in this proportion was the same in Texas as in California?


title: Does relative Kullback-Leibler divergence exist?

text: Suppose I have two multivariate normal distributions. I have computed the KL divergence ($d_{LK}(N_1, N_2)$). Is there a way to measure a relative divergence between these two distributions? For instance in a deterministic case, there is absolute error and relative error. Let's say if KL divergence is analogous to absolute error, is there an analogous case for relative error. Any references to literature would be highly appreciated!


title: Can I regress an index value with variables used to create the index?

text: I developed an index value (vulnerability score scale of 0 to 1) using a series of variables. I would like to regress these variables with the index value to determine the relative predictive power of each variable. Can I do this? I ran the regression and came up with standardized B coefficients. I then interpreted those as relative contribution of the variable towards predicting the index value (vulnerability score). I know one cannot regress a variable against itself but I am essentially doing this but primarily am just looking at determining to what degree each variable predicts the indexed value. Any insights would be helpful if this is a proper use of regression or if there is an alternative method to assess this. Thanks!


title: Statistics Proof that $E[(X-Y)^2] = 0$

text: If X and Y are standardized variables and are perfectly positively correlated with respect to each other, how can i prove that $E[(X-Y)^2] = 0$?


title: Using weighted regression to obtain fit lines for which I only have summary data

text: I'm working with a dataset for which I only have means, standard deviations, and sample sizes for different levels of a continuous predictor. E.G.               Y  X SD_Y N_Y     5  1   3   4     10 2   6   2     15 3   2   8      I would like to determine the regression line that fits this data. I'm wracking my brains to remember how data points should be weighted in a linear regression (I'm also interested in using a generalized linear model as well)- by sample size, variance, SD? Any pointers?


title: Regression Proof that the point of averages (x,y) lies on the estimated regression line

text: How do you show that the point of averages (x,y) lies on the estimated regression line?


title: Issues in estimation and plot

text: I am learning adaptive filters and testing the performance of using Least Squares and Kalman filter for parameter estimation for $y = X + \text{noise}$. The model is autoregressive AR(2) model $$y(t) = ay_{t-1} + by_{t-2}$$ and a Moving Average MA(2) model. I have obtained the parameters $(\hat{a},\hat{b})$ for AR and MA. I want to show the plot of the desired signal vs the actual signal so that it shows how close the estimated signal is from the actual. How do I recover the signal and show the plot?


can everyone tell me how to use ggplot2 make recursive partitioning survival tree? I know " party " package can make a plot in R, however the plot is not looking good, and all the plots in the rest of my report are made by ggplot2. for my plot, I printsreen and added labels manually. if there is a way to input data and plot in R also works, Thanks!! ![enter image description here](http://i.stack.imgur.com/pcbhx.png)


title: Is there any MOOCS on statistical analysis of Applied behavioral economics?

text: I couldn't find a single MOOCs on Applied behavioral economics and statistics, on how can we used statistical analysis in behavioral economics. If someone can point me in the right direction, it would be great


title: Assessing the accuracy of zero-inflated beta regression models

text: I have fitted a zero-inflated beta regression model to my data in R, using the gamlss package. However, I am unsure of how to assess the fit of the model to my data, i.e. finding a coefficient of determination. Does anyone know if this can be easily calculated from a zero-inflated beta regression model? Thanks in advance


title: Investigating relationships between variables- multiple and simple linear regression

text: please help a sociologist struggling to get to grips with R and statistics in general..! I've got a data set with 11 variables. I need to investigate the possible relationships between one of the variables (percentage of smokers in a population) and the rest of them- things like unemployment rate, education level, and so on. I'm working purely with linear regression. I'm getting there with correlation, but having looked around at numerous examples on the web, I'm still unsure with regards to method. Apologies if it's completely obvious, but I can't seem to find an answer anywhere. As I see it, it seems sensible to look for correlation between smoking and all the other variables firstly with scatterplots, and then by calculating the correlation coefficient if there is an identifiable linear relationship (Pearson's r or Spearman, depending on normality). If we then continue and do separate simple linear regressions between smoking rate and our identified variable, that all seems well and good. But where does multiple regression fit in here? Does it make any sense, statistically speaking, to do simple linear regression first and then multiple linear regression also? Or would it be best to just jump straight from testing correlation to multiple regression? Any help much appreciated!