title: Unbiased estimate of population standard deviation: is sqrt(2) a superior correction?

text: # Background on bias correction constants The standard deviation is calculated like this: $$ SD = \left(\frac{1}{N-constant} \sum_{i=1}^N (x_i - \overline{x})^2\right)^{1/2} $$ Following Wikipedia's entry on the standard deviation, the biasedness of the estimation of a population SD given a sample depends on $constant$ in the following way:   * $correction=0$: just sample sd, heavily biased to be smaller than population sd.   * $correction=1$: Bessel's correction, less biased but still smaller.   * $correction=1.5$: "rule of thumb", the best single value for an unbiased estimate. # Question I simulated this and found that $\sqrt{2}$ is an even better value for $constant$, especially for small samples (n < 10) where $1.5$ overestimates population SD. Have I just discovered something fantastic or am I missing something here? # Simulation For each of the sample sizes $n=2,3,5,10,20,60,100,200$, I generated 3.000 samples using `rnorm(n, 0, 15)`. For each sample size, I then estimated population SD using each of the above constants. Here's the result: ![simulation results](http://i.stack.imgur.com/KdLOl.png) Each plot is a different estimation constant. The "error" in the title is `mean(sd.estimations - sd.real)`. The red line is the real SD. The blue line shows the estimated sd. The vertical gray lines mark change in sample size. Points show individual sd-estimations. It is clear that $\sqrt{2}$ is superior to $1.5$. This is true for large sample sizes as well, even though it's not clear from this plot. Here's the R script that generated these plots. # Update and conclusion $\sqrt{2}$ is close to the analytically correct solution but doesn't outperform it. It remains a heuristic which could be used out of laziness or for computational efficiency with small sample sizes. Actually, the closest approximation depends on the sample size you want to calculate for. Here's a few optimal values for different sample sizes:   * $4 < n < 10$: $\sqrt{2.15} = 1.47$ deviates with 0.4% at most.   * $10 < n < 50$: $\sqrt{2.22} = 1.49$ deviates with 0.04 % at most.   * $40 < n < 300$: $\sqrt{2.2465} = 1.4988$ deviates with 0.0025% at most. As sample size increases, the constant approaches the "1.5 rule of thumb". Thus the conclusion is that $\sqrt(2)$ is quick-and-dirty for small sample sizes. For larger samples, reasonable approximations can be made with 1.5. And just to be clear: Bessel's correction is still the right way to get unbiased when estimating _variance_. The observations above only pertain to the estimation of the _population standard deviation_.


title: Regression forest: Leaf node and information gain

text: Regression forests are basically random forests, however used for regression. They basically use the same framework as decision forests use for classification with a few parts exchanged. Two of these parts are:   1. How to calculate the information gain when deciding how to split at a node?   2. How do we make a prediction when we reach a leaf? An interesting report (however 150 pages strong) concerning this topic is the following: http://research.microsoft.com/apps/pubs/default.aspx?id=155552 This report touches the above mentioned problems, I however do not really understand the possible solutions. 1) Page 52, Appendix A. They define the continuous formulation of information gain to be: $ I_j = \sum_{v \in S_j} log(|\Lambda_y(v)|) - \sum_{i \in \\{L,R \\}} (\sum_{v \in S_j^i} log(|\Lambda_y(v)|)) $. As I dont really know where else to ask: Does anyone know what || is in this case? How would you implement a function to compute I? 2) When we reach a leaf we have a number of data points reaching this leaf, say of the form (x,y). From my understanding we would now fit a function f(x) to these data points (e.g. by least squares) and would then use f(x) to predict an y for a future testing point x. In the report (page 50) they talk about a probabilistic linear fitting, that returns a conditional distribution p(y|x). How could that be done? I dont understand how we can compute a probability density function here. I hope anyone is familiar with the above mentioned report and might be able to help :) Thank you!


title: Feature importance scores of SVM multiclass one-vs-one design

text: Info about dataset: 5 classes, 200 trials, 100 features. (I know about the trial to feature ratio being very low, but can not avoid this here and still got well enough classification results.) Within a ten-fold cross-validation a SVM classifier was trained for each fold (linear nu-svm, optimal nu was identified by inner ten-fold cross-validation). To handle the multiclass problem the libsvm one-vs-one approach is used: $10=5*4/2$ binary classifiers are trained and then a majority voting scheme used. If two or more classes get the same maximum number of votes one of these is choosen at random (at this point I modified the libsvm default behaviour). * * * I wanted to compute feature importance scores. I tryed the following procedure:   * load the model of one fold and the matching test data   * predict the test data and remeber the prediction accuracy as baseline   * repeat the following for every feature 500 times:      * permute one feature and get the prediction accuracy on the altered data   * In the end: for every feature compare the prediction accuracy for the intact data with the ones based on the altered data and calculate the average decrease in accuracy. Averaging these values over all folds should give me an overall importance score i.e. an overall average decrease in accuracy for every feature. The higher this value the more important the feature. * * * I thought this idea should work out, but the results of this do not look sensible to me: For more or less every feature this yields an average _increase_ in accuracy when permuting this very feature. At least for me that does not make really much sense. Further more none of these averaged values seems to be significant in the sense, that there is no feature for which more than 95% percent of the accuracies with this feature permuted are lower than the baseline. Any suggestions how to get some sort of importance scores or how to improve my method? (In Matlab this is actually really really slow procedure as I implemented it at the moment. :-()


title: How to include GAM response error in a Monte Carlo simulation

text: I am running a Monte Carlo simulation, using the results of a GAM (response) as the basis for my overall model. I would like to incorporate the error in the GAM into the final result. Since predict.gam (mgcv) includes se.fit, is there a way to resample the response around this standard error? I was thinking of using a random normal distribution based on the mean response and the standard error, but I don't know if that is valid, given that the standard error is calculated from a covariance matrix and not a population (code below). The R2 of the GAM is something like .87; is there a standard way to account for this unaccounted for variance in the response? My pseudo code (R):               library(mgcv)      model = gam(B ~  s(C) + s(log(D)), data=data.B, family=gaussian,                  link="identity")       loop:               a= predict.gam(model, newdata = input, type="response",                    na.action=na.pass, se.fit=T)           a_se = rnorm(a$fit, mean=a$fit, sd=(a$se.fit*(????)))       Many thanks


title: Do you reject the null hypothesis when $p < \alpha$ or $p \leq \alpha$?

text: This is clearly just a matter of definition or convention, and of almost no practical importance. If $\alpha$ is set to its traditional value of 0.05, is a $p$ value of 0.0500000000000... considered to be statistically significant or not? Is the rule to define statistical significance usually considered to be $p < \alpha$ or $p \leq \alpha$??


title: SVM confidence according to distance from hyperline

text: For a probabilistic multi-class classifier we can get probabilities of membership of a new point $x$ to each class $y_i$; in case of 3 classes suppose that we get $P(y_a|x) > P(y_b|x) > P(y_c|x)$, thus the most probable class of x is $y_a$. Now suppose that we have a multi-class svm where we can get scores of membership of $x$ to each class (according to distances from hyperlines); in case of 3 classes suppose that we get $Score(y_a|x), Score(y_b|x), Score(y_c|x)$, How is in this case the first, second and third most likely class of $x$ (without converting these scores to probabilities) ? Usually I get positive and negative values like for instance $Score1 = -8622, Score2 = 5233, Score3 = -665$


title: Taking average of linear regression results

text: Say I have a very large sample of observations, like a thousand of them, and I want to do (multiple) linear regression on them, i.e., regressing the response variable on the predictors. Now instead of doing so, I would randomly sample a subset of that large sample (say half of it, or a third of it, or even a fourth) and do regression analysis on that subset only. So I do this for $n$ times and each time I have a set of results, like coefficient estimates, $p$-values, etc. Finally I take the average of all the $n$ independent sets of results to get the final regression results. The question is: Would this rather "awkward" method be equivalent to the straightforward one (doing regression on all the observations that I have)? By "equivalent" I mean the coefficient estimates (signs and magnitudes) and $p$-values resulted from the two methods are the same, or insignificantly different. My second question is, as we all know, that in order to have a good regression analysis, it is often required to have a sufficiently large number of observations. What about the other way round? If I have too large a sample, like thousands of observations, should I just fit them all or should I sample a subset from it and fit?


title: Solve a statistical equation

text: I got the following exercise where i have to understand the steps theres are done: Since $z_T \sim N(0,T)$: $$\text{Prob}((\mu-0.5\sigma^2)T+\sigma z_T>rT)$$ $$\text{Prob}(z_T>-\dfrac{(\mu-r-0.5\sigma^2)T}{\sigma})$$ $$\text{Prob}(z_T<\dfrac{(\mu-r-0.5\sigma^2)T}{\sigma})$$ Then i don't understand how i get to this: $$\text{N}(\dfrac{(\mu-r-0.5\sigma^2)T}{\sigma\cdot\sqrt{T}})$$ $$\text{N}(\dfrac{(\mu-r-0.5\sigma^2)\sqrt{T}}{\sigma})$$ I hope someone can help - and thanks for your time!


title: Converting R's Random Forest model to Mahout

text: I need to implement an RF created by R on Hadoop. I know a way is to recreate the random forest with the same training data using Mahout. But is it possible to convert the RF created by R directly into a format that Mahout can use, for scoring?


title: Regression in SEM programs vs regression in statistical packages such as SPSS

text: **Contex** Multiple regression with continuous variables. Conventional statistical packages, e.g. SPSS, lm() in R, typically give me an _F_ value, dfs, and a significance test on whether the model is performing well. I can also use ANOVA to compare whether one model is significantly better than another by finding out whether the increase in R square is significant. Structural equation modeling programs, e.g. Mplus, lavaan package in R, also allow me to run regression with added benefits of missing data handling through full information maximum likelihood. Even though parameter estimates and R squares are similar between these two approaches, I do not get the _F_ value and the significance testing of model doing regression in SEM style. Instead, I get model fit information with df being `0`, and I cannot figure out how I could test whether one regression model is better than another. **Questions**   1. What are the differences between these two approaches, if any?    2. In SEM style, how do I know a regression model is good, and how can I compare whether one regression model has significantly larger R square than another?


title: Detection of periodic patterns using DWT

text: Is it possible to detect a periodic pattern in a time series using discrete Wavelet Transform? Is there any package in R to do this job?


title: How can I estimate a score for my users, on if they will match ad campaign criteria?

text: I have ad campaigns that I want to display to users. Each one has criteria a user must meet in order for them to be valid for the campaign. For example, "User must be in the USA", or "User must like the albums by at least one of the artists [The Beatles, Michael Jackson], AND user must be in USA". Right now I just let the user see all of the ads as they browse though my site, but I want to optimize for certain scenarios where I have scarcity (only show one ad on the home page, or one featured ad in search results). I want to be able to pick the best ad to display. At least half the ads have such custom requirements that I can't apply the criteria to see if they are valid to target. After the user clicks the ad placement, they are asked to fill out more data, at that point we have enough info to know if they are valid, and we match the user with the advertiser or product. We then collect revenue. I am a software developer, and I don't know much about stats or machine learning. I want to know if there is some kind of machine learning or statistical model that I can use to help optimize this process. I'm imagining some kind of system where I run our ad targeting criteria against our entire user base's data (we have over 5 million users), and generate a set of tuples:               {user: John Doe, ad_campaign: Coca-Cola, valid_to_target: Yes/No/Maybe(but not enough info)}      Then, I would feed that set of tuples as a training set into some kind of ML model (not sure if I'm using the right terminology). That process would result in a model, a piece of code I can execute in my ad service, that would be able to take as input a brand new user, and output a set of tuples similar to the training data tuple above, except these tuples would also have a score indicating how likely that user is to be valid/targetable. I guess a `Yes/No/Maybe` could also work, but sometimes the criteria involve several user data properties, so I would want the model to be able to give a higher score for users who look like they might be able to be valid, maybe have data for 5 fields out of 6 that are part of the ad targeting criteria, and the 5 fields they have data for don't rule them out for that campaign. So that user would see that ad over one where they have maybe 1 field out of 4. Then, I could apply the ad revenue/costs with the scores to know know exactly what to show the user. I'd like to know if this is possible, and if so, what are some canonical references that would help me learn to accomplish it. I see that there are libraries like `Infer.NET` (I'm a .NET developer mostly) that might be able to help.


title: Choosing clusters for k-means: the 1 cluster case

text: Does anyone know a good method to determine if clustering using kmeans is even appropriate? That is, what if your sample is actually homogenous? I know something like a mixture model (via mclust in R) will provide fit statistics for the 1:k cluster case, but it seems like all of the techniques to evaluate kmeans requires at least 2 clusters. Does anyone know of a technique to compare the 1 and 2 cluster cases for kmeans?


title: Rounding when making a histogram

text: Suppose that when making a histogram, one encounters a datum at a bin boundary. Is there a convention on how to round it? For example, suppose my data are integer percentages, running from 0% to 100%. I want to make a histogram with boundaries at 10%, 20%, etc. Of course I could label the histogram in a way to avoid the issue, e.g. by labeling the categories 0-9, 10-19, etc. But suppose I don't want to label them like that. I could save space on the x axis by simply labeling the bin boundaries of 10, 20, 30, etc. But then if I encounter a datum at 20, is there a convention as to whether it goes into the 10-20 bin or the 20-30 bin?


title: What type of time series model would be good?

text: I am trying to fit a time series model to the following data. It seems to be seasonal. Would an ARIMA model be good? ![enter image description here](http://i.stack.imgur.com/7LaZu.jpg) Here is the data:               Count      2 1 4 5 4 8 7 11 4 4 11 7 10 7 0 19 13 13 11 9 8 16 10 12 9 7 21 9 10 6 7 19 18 9 19 15 14 17 9 10 10 13 15 20 15 12 15 16 The numbers are separated by spaces.


title: Confusion related to log transformation

text: I found this preprocessor where they had the data log transformed with $log_{10}(100+data)$. I am not sure what this log transformation exactly does. Any suggestions?


title: finding global minima of a random forest estimator

text: I have a random forest regression model with 1000 trees, having 16 parameters (using python scikit-learn). The estimator can predict a target value with cross validated r2 score of 0.87 +/- 0.03. I would like to find the global (or at least close to global) maximum of the random forest, in order to identify the best values of the 16 parameters that yield this optimum. What is the best algorithm to reach the maximum of a random forest estimator? Is there a python implementation of this algorithm?


title: Negative variance from inverse Hessian matrix

text: I used nlm function in R to do the optimization. When I calculated the correlation between estimated parameters using the inverse of Hessian matrix, I got negative values on the diagonal. My questions are what could be the reason for that issue, and what could be possible solution?


title: Can LDA assign more than one topic for a word?

text: I have just started reading about Latent Dirichlet Allocation LDA and want to apply it to my project. May I know if LDA is able to assign a topic to more than one word? For example, Article A talks about "river banks" while Article B talks about "The role of banks in finance". Hence, will LDA allow the word "banks" to potentially be assigned to two different topics?


title: US heatmap using zipcodes

text: We are attempting to create a US heatmap using only zipcode and corresponding values using R. While there are a number of examples that use lat and long data for creating such maps, is there one that will accept zipcodes and colour the corresponding regions on a map ? Preferably using some package that will produce presentation quality images (ggplot2, etc). Thanks in advance, - Raj.


title: Scatter plot or 2D histogram for mixture of Gaussian fitting?

text: I am ust a little bit confused about whether I should use a scatter plot or a 2D histogram to fit a mixture of Gaussians. The problem is I have intensities from one image with 2 channels. For example, a RGB image but using just Red and Green channels. Then I plot the scatter plot and try to do a mixture of Gaussian fitting for the data points in the scatter plot. However, thinking about the problem again, I am not sure if I should do a Mixture of Gaussian for the 2D histogram. Essentially, a 2D histogram is different from a 2D scatter plot, am I right? Sorry about my confusion and thanks for your help.


title: Statistical significance between groups of a summation variable

text: The general problem is that I have an calculated variable that represents agreement across parts. The variable is a summation in the form (part/count of parts)^2 for possible parts. The summation variable is between 0 and 1. How would I go about testing for significant differences between calculated part agreement for two study conditions? The specific situation involves gestures that study participants used to accomplish a particular task. Here, agreement means the number of times a particular gesture was used over the total number of gestures for each task (squared as per above)--summed for all gestures per task. I would like to calculate whether or not agreement differs significantly per gesture across two groups (study compares two conditions) and am unsure how to approach this problem.


title: Given a series of past coin tosses from an unfair coin, how can I calculate the confidence for the next toss result?

text: So this is really difficult for me, but I would like to know if this is possible. Let's assume I have an unfair coin (absolutely no assumptions can be made about the coin), and my past result has been `T-T-T-H-T-H-T-T-T` and I want to predict the the result of the next toss. In a way that let me formulate my answer in this way: "There is an X% confidence that the result of the next toss will be T". So I am trying to solve for X.


title: Whether distributions with the same moments are identical

text: Following are similar to but different from previous posts here and here   1. Given two distributions which admit moments of all orders, if all the moments of two distributions are the same, then are they identical distributions a.e.?   2. Given two distributions which admit moment generating functions, if they have the same moments, are their moment generating functions the same? Thanks!


title: Calculate Similarities between two arrays with different scales

text: In creating a content-based recommender, I have several binary measurements, combined with integer measurement values and I wish to compare the two vectors. For example, consider the following two users and traits t1 - t5, where the traits could be, for example:   * meets_requirement1   * meets_requirement2   * meets_requirement3   * years_experience,    * years_college `... [ t1 t2 t3 t4 t5]`               User1   0    1    1   18    4     User2   1    1    1   15    2     User2   1    1    0    4    8      The data are (clearly) not Gaussian across rows. How do I compute similarity scores between these users? Do I need to normalize down traits first? I'm looking for a result like the following (in Python, just to show the result. I don't need it in Python.):               a = array([[0, 1, 1, 18, 4],     [1, 1, 1, 15, 2],     [1, 1, 0, 4, 8]])          print squareform(1 - pdist(a, metric='cosine'))      print squareform(1 - pdist(a, metric='correlation'))     print squareform(1 / (1 +  pdist(a, metric='euclidean')))      Yielding the output:               [[ 0.     0.994  0.627]      [ 0.994  0.     0.566]      [ 0.627  0.566  0.   ]]          [[ 0.     0.991  0.384]      [ 0.991  0.     0.273]      [ 0.384  0.273  0.   ]]          [[ 0.     0.211  0.064]      [ 0.211  0.     0.074]      [ 0.064  0.074  0.   ]]      But I want to be sure I'm not over emphasizing the t4 & t5 values.


title: Clustering spatio-temporal data?

text: I have data in the form of timestamp,lat,long which is gps data for users. I'm new to data mining and want to understand how can I start clustering these data to understand more about it. Should I like build a matrix of one trajectory v/s the other based on some distance metric and then apply some clustering algorithm on it? Data will be a trajectory for each user. Each user will have a sequence of points in the form of (timestamp,lat,long) starting at point A to point B. I want to cluster the trajectories.


title: combining and contrasting time course GLMs using R

text: I am analyzing some time course data in which I have set up a GLM using R for each subject. Each GLM I want to run is an attempt to extract estimates of different behavioral conditions effects on the time series, and each GLM is all within a single individual.               glm(dilation~happy+fear+neutral+post_decision, data=dilation_data)      the matrix looks similar to this, only it has about 50,000 time points (so the same number of rows):               dilation = intercept happy fear neutral post_decision             0.3          1        0     0     0        0           0.5          1        0.2   0     0        0         0.6          1        0.8   0     0        0          0.9          1        1     0     0        0           0.5          1        0.4   0     0        0.2         0.8          1        0.2   0     0        0.4          1.2          1        0     0     0        0.8          0.9          1        0     0     0        1             1. Is it going to be a problem in terms of contrasting different estimates or the model in general if my time courses of the regressors slightly overlap (see happy and post_decision in example above), do they need to be made orthogonal for the glm to be correct? If so what is the best way to do this using R.   2. I have 40 subjects, my plan is to run a glm for each subject and combine and contrast the different outputs. What is the best way in which I can contrast the difference between the GLM outputs for happy-neutral for example, and is it recommended to do this before or after combining subjects?   3. What is the ideal way for me to combine the 40 GLMs, so I have an average estimate of the change from baseline for happy lets say, and what test would be ideal for trying to say happy is significantly different than neutral based on this initial GLM for all 40 subjects?


title: Generate random number from a piecewise exponential distribution

text: I would like to generate a random number from a piecewise exponential distribution. I consider that the time-scale is divided in $J$ intervals with bounds $(s_{j-1},s_j]$, for $j=1,...,J$, and corresponding rates $\lambda_j$. Considering the memoryless property of the standard exponential distribution, is it correct to generate the random number with the following procedure for $i=1,...,J$?   1. Generate $x_i$ from a standard exponential distribution with rate $\lambda_i$;    2. If $x_i<s_i$ the random number is obtained as $y_i<x_i+s_{i-1}$, otherwise continue with the next $i$.  In addition, the PDF of this piecewise exponential distribution is given by: $$ k(t)=\prod_{h=1}^{j-1}(e^{-\lambda_h(s_h- s_{h-1})})(\lambda_j)(e^{-\lambda_j(t-s_{j-1})})I(s_{j-1}<t\leq s_j) $$


title: Comparing multiple groups to a control group

text: I've made several biochemistry experiments, and in which I compare several groups to a control group. I want to compare the results of 8 such experiments, so basically I have 8 measurements in each group, but the measurements in the control group is all 100%. How do I compare them? I.e., I have four groups, the first three are 75%, 45.5%, 50% etc., but the control group to which I compare them is all 100%


title: How to calculate post-hoc power analysis for multiple regression models?

text: I have few questions that came up while I was trying to calculate post hoc power analysis for a regression model:   1. I calculated post hoc power analysis for a regression model that produced a statistically significant (adjusted) R2 = .36 [F (2, 433) = 121.2, p < .001]. N=405. Number of predictors=2. I used F-test, statistical multiple regression (fixed model, R2 deviation from 0) and received a very low power - 0.05. Does that make any sense?    2. I have analyzed data by performing a regression model in three steps - (a) two psychological independent variables, (b) adding a political independent variable (3 variables), (c) adding a political demographic independent variables (five variables). Am I supposed to report power analysis for each of the steps? Or - should use F-test, statistical multiple regression (fixed model, R2 increase)?


title: Retrieving standard deviations from OLS estimates

text: I am running a fixed effects regression in Stata: $y=\alpha+\beta D$ (omitting the FE), where $D$ is a dummy variable. Basically, I want the mean of group $D=0$ and of group $D=1$ after controlling for time invariant heterogeneity via fixed effects, along with their standard deviations. Stata reports coefficients and standard errors. Now, the means after FE are simply $E[Y|D=0]=\alpha$ and $E[Y|D=1]=\alpha+\beta$, because the intercept in Stata fixed effects regression is just the mean. But how do I get the standard deviations? Groups may not have the same number of observations. If I had just one group, I would estimate $Y=\alpha$, and from the standard error of the intercept I can calculate the sd with $SD=SE*\sqrt{n}$, where $n$ is the number of observations. But that doesn't seem to work with two groups and the dummy. I tested with a standard OLS (without fixed effects) and, say, the SD of group $D=0$ cannot be computed via $SD_0=SE(\alpha)\sqrt{n_0}$ (a degrees of freedom correction of 2 doesn't change much). So how can I do it, get the SD of "FE corrected" means?


title: Evaluating if the new process was better than the previous?

text: Initially things were done in a particular way (A) and I changed to in some key ways to do things a bit differently (B). Now I wish to evaluate whether the change was for the better. What tests can I do to validate this claim? I have historical and observational data for 3 years on each of the approaches. Which statistical tests can I use to measure/quantify the effectiveness of the change? For example: Items were selected for implementation were selected in a particular manner. I changed the selection process to be more rigorous and value-centric (for some notion of value) and claimed that the items selected by my approach are superior with respect to their value for the organization/implementation. I now wish to conduct a retrospective analysis to see if things are better as a result of my approach or just due to chance. I was suggested to use MANCOVA which I'd never heard of/used before. I'm working on understanding it but is there something else that I should look at? How are process changes statistically evaluated 'in the wild'? (Or just left to observational data?)


title: Should you use normalized or non-normalized data to develope your model?

text: I am developing a linear model with 13 variables, including the target variable (online purchase revenue for items). So, I first built model1 with regular variable and then build model2 after normalizing of the data. I have copied the coefficients for two models here :               Model1(Without Normalized Data)     Coefficients:                   (Intercept)                        xid                   xcartadd                       6.386e+01                 -4.301e-03                 -1.229e+02                    xcartuniqadd            xcartaddtotalrs                xcartremove                       1.239e+02                  7.788e-02                 -1.424e+02               xcardtremovetotal        xcardtremovetotalrs              xproductviews                       5.588e+02                 -3.445e-02                  1.369e+01                   xuniqprodview              xprodviewinrs       xsizeselecteduniview                      -1.530e+01                  5.401e-04                 -1.299e+02         xsizeselectedtotalviews  xsizeselectedtotalviewsrs                       6.280e+01                 -2.453e-02            Model 2(With Normalized data)     Coefficients:                     (Intercept)                          xid                         3.900e+02                   -4.301e-03                        xcartadd_n               xcartuniqadd_n                        -2.623e+03                    2.069e+03                 xcartaddtotalrs_n                xcartremove_n                         1.785e+03                   -1.721e+02               xcardtremovetotal_n        xcardtremovetotalrs_n                         4.474e+02                   -5.360e+01                   xproductviews_n              xuniqprodview_n                         7.979e+03                   -7.378e+03                   xprodviewinrs_n       xsizeselecteduniview_n                         4.757e+02                   -1.218e+03         xsizeselectedtotalviews_n  xsizeselectedtotalviewsrs_n                         1.044e+03                   -5.374e+02      **and my questions are :**   1. Is it appropriate to take only normalized data into model or non normalized data?   2. Is it appropriate to take combination of normalized data as well as non normalized data?   3. How can I choose most appropriate predictor variable from them for model?


title: Fit a sinusoidal term to data

text: Although I read this post, I still have no idea how to apply this to my own data and hope that someone can help me out. I have the following data:               y <- c(11.622967, 12.006081, 11.760928, 12.246830, 12.052126, 12.346154, 12.039262, 12.362163, 12.009269, 11.260743, 10.950483, 10.522091,  9.346292,  7.014578,  6.981853,  7.197708,  7.035624,  6.785289, 7.134426,  8.338514,  8.723832, 10.276473, 10.602792, 11.031908, 11.364901, 11.687638, 11.947783, 12.228909, 11.918379, 12.343574, 12.046851, 12.316508, 12.147746, 12.136446, 11.744371,  8.317413, 8.790837, 10.139807,  7.019035,  7.541484,  7.199672,  9.090377,  7.532161,  8.156842,  9.329572, 9.991522, 10.036448, 10.797905)     t <- 18:65      And now I simply want to fit a sine wave $$ y(t)=A\cdot sin(\omega t+\phi) +C. $$ with the four unknowns $A$, $\omega$, $\phi$ and $C$ to it. The rest of my code looks is the following               res <- nls(y ~ A*sin(omega*t+phi)+C, data=data.frame(t,y), start=list(A=1,omega=1,phi=1,C=1))     co <- coef(res)          fit <- function(x, a, b, c, d) {a*sin(b*x+c)+d}          # Plot result     plot(x=t, y=y)     curve(fit(x, a=co["A"], b=co["omega"], c=co["phi"], d=co["C"]), add=TRUE ,lwd=2, col="steelblue")      But the result is really poor. ![Sine fit](http://i.stack.imgur.com/IS0ae.png) I would very much appreciate any help. Cheers.


title: RTextTools for coded sentences

text: I'm new to text mining and I'm not sure if could be applied here: I have labeled sentences but the words I work with are codes, so each word length ranges from 5 to 15 letters and there are only 4 different letters (A, B, C, D). Looks like this :               Sentence1: AAABCD BAAAA ABCCBAA CCCCC...... AAAAAAABB CCCCCCCAAA (Positive labeled) ...      Sentence10000: ACCCCCCCCCCCD ... AACCDD (Negative labeled)      Each sentence could be composed of a different number of words in a different order, ranging from 6 to 50 words. Could RTextTools work in this case? When the matrix is created you have to select a language and I have no language.


title: post-hoc test on mixed-effects model

text: I have an experiment with bacteria for which I measured growth curves of several bacterial clones in two different types of media (treat). I then extracted several growth parameters from the curves and am analysing one of these: lag phase (lambda). clone has 13 levels and treat has 2 levels. I now want to know, for each clone, if my lambda is different between the two treatments. As I carried multiple non-independent replicates for each clone, I created a mixed effects model of the following form               m1<-lme(lambda~treat*clone,random=~1|replicate)      As the interaction is significant, I used lsmeans to try to get the comparisons between the two treatments for the different clones. Below is the code I used for the function and the output. Now there are two things puzzling me: 1) why is lsmeans unable to estimate the df?; 2) Why is the SE the same across most of the different clones?               lsmeans(m1,pairwise~treat|clone,adjust="tukey")                    $`treat:clone lsmeans`        treat clone    lsmean       SE df asymp.LCL asymp.UCL          glu     1  274.5716 40.05260 NA  196.0700  353.0733      glu.ALA     1  316.6030 40.05260 NA  238.1014  395.1047          glu    11  628.6852 40.05260 NA  550.1835  707.1868      glu.ALA    11  301.2156 40.05260 NA  222.7139  379.7172          glu    19  586.9249 40.05260 NA  508.4232  665.4265      glu.ALA    19  622.3766 40.05260 NA  543.8749  700.8782          glu    26  712.6377 40.05260 NA  634.1360  791.1393      glu.ALA    26  350.8608 40.05260 NA  272.3592  429.3625          glu    28  491.2373 40.05260 NA  412.7357  569.7390      glu.ALA    28  314.1701 40.05260 NA  235.6685  392.6718          glu    31  559.2341 40.05260 NA  480.7325  637.7358      glu.ALA    31  452.9919 40.05260 NA  374.4902  531.4935          glu    40  569.0526 40.05260 NA  490.5509  647.5542      glu.ALA    40  592.2210 40.05260 NA  513.7194  670.7227          glu    41  430.7679 40.05260 NA  352.2663  509.2696      glu.ALA    41  456.9167 40.05260 NA  378.4150  535.4183          glu    48  739.8913 40.05260 NA  661.3896  818.3929      glu.ALA    48  709.2002 40.05260 NA  630.6985  787.7018          glu    51  388.9805 40.05260 NA  310.4788  467.4821      glu.ALA    51  395.3881 40.05260 NA  316.8865  473.8898          glu    52  499.6519 40.05260 NA  421.1502  578.1535      glu.ALA    52  471.2882 43.87542 NA  385.2940  557.2824          glu    61 1200.5866 40.05260 NA 1122.0850 1279.0883      glu.ALA    61  817.7402 40.05260 NA  739.2385  896.2418          glu   Anc  262.7183 25.33149 NA  213.0695  312.3671      glu.ALA   Anc  302.1657 26.22058 NA  250.7743  353.5571          $`treat:clone pairwise differences`                           estimate       SE df  z.ratio p.value     glu - glu.ALA | 1   -42.031407 56.64292 NA -0.74204 0.45808     glu - glu.ALA | 11  327.469629 56.64292 NA  5.78130 0.00000     glu - glu.ALA | 19  -35.451732 56.64292 NA -0.62588 0.53141     glu - glu.ALA | 26  361.776858 56.64292 NA  6.38697 0.00000     glu - glu.ALA | 28  177.067207 56.64292 NA  3.12603 0.00178     glu - glu.ALA | 31  106.242222 56.64292 NA  1.87565 0.06073     glu - glu.ALA | 40  -23.168464 56.64292 NA -0.40903 0.68253     glu - glu.ALA | 41  -26.148743 56.64292 NA -0.46164 0.64435     glu - glu.ALA | 48   30.691103 56.64292 NA  0.54183 0.58795     glu - glu.ALA | 51   -6.407644 56.64292 NA -0.11312 0.90994     glu - glu.ALA | 52   28.363655 59.40760 NA  0.47744 0.63306     glu - glu.ALA | 61  382.846452 56.64292 NA  6.75895 0.00000     glu - glu.ALA | Anc -39.447421 36.45824 NA -1.08199 0.27928         p values are adjusted using the tukey method for 2 means


title: Standard error of probability estimate

text: Take a look at Table 1 on page 268 in http://www.math.ku.dk/~rolf/teaching/thesis/DixonColes.pdf It says at the end of the previous page that the standard errors are computed "on the basis of an underlying multinomial model". The sample size is 6629 (also mentioned on the previous page). If I understand things correctly the standard error should be $\sqrt{\frac{p\cdot(1-p)}{n}}$, but computing them in this way I can't have them match the results in the paper. For example the entry 33.4 (0.74). I get that the standard error should be $100 \cdot \sqrt{\frac{0.334 \cdot 0.666}{6629}} = 0.57927707090437$ (the factor 100 just to norm to percentages). Have I misunderstood something here?


title: Can I compare models from linear regression and nonlinear regression using RMSE?

text: I am comparing multiple published equation forms, refit with independent data. I'm trying to be true to the original authors' methods as much as possible. Therefore, I have 3 linear equations (fit in R using lm()), two of which use transformed Y-variables, and one equation fit using nonlinear regression (fit in R using the gnls() function). In all instances cases I'm weighting the residual variance structure using the inverse of one of the predictors to account for observed heteroskedasticity. I have been evaluating the models using R2, and RMSE- using back-transformed data for the two models with transformations. I've calculated RMSE "by hand" using the following equation:                RMSE<-sqrt(sum(residuals(Equation)^2)/length(residuals(Equation))-2))      Should I use similar code to calculate RMSE for the linear and nonlinear regression models? Is the metric still a valid statistic for comparison, or am I missing some important assumption? Edited: I initially stated that I was also comparing models using AIC; I later recalled that AIC would not be appropriate if the Y-variables were transformed because the models would be estimating different things.


title: How to interpret "weight-position" plot when using self-organizing map for clustering?

text: I used MATLAB neural network toolbox to train a self-organizing map for a given data set. The obtained "weight-position" plot is given as follows. I do not think this plot looks good in comparison to the sample plot given in the MathWorks website. ![enter image description here](http://i.stack.imgur.com/AlZEr.png) I followed the standard MATLAB routine to train this SOM. Thus, I am having three questions on this case study.   1. Which kind of information can I obtain from this plot?   2. Given the shape of this plot, what can I say for the input?   3. Since the training method is trivial and standard, I am wondering whether SOM can fit to my data set? My data set is composed of 3,000 data points, and each data points has 20,000 dimensions. The example given by MATLAB is only for 4 dimensions.


title: Invertibility of AR(p) model

text: Notation: $\dot{Z}_t = Z_t - E(Z_t)$, so that it is centered at 0. $a_t$ stands for the residual and we assume the $a_t$ are independent and normally distributed with mean 0 and constant standard deviation $\sigma_a^2$. And, when I say stationarity, I mean weak stationarity and not strict stationarity. In my time series class, we have been given that an MA(q) model is of the form $$\dot{Z}_t = \theta_q(B) a_t$$ where $$\theta_q(B) = 1 - \theta_1 B - \theta_2 B^2 - \cdots - \theta_q B^q.$$ We were told MA(q) is always stationary and it is invertible when all of the roots of $\theta_q(B)$ are outside the unit circle. For an AR(p) model, we were told it is of the form $$\phi_p(B) \dot{Z}_t = a_t$$ where $$\phi_p(B) = 1 - \phi_1 B - \phi_2 B^2 - \cdots - \phi_p B^p.$$ We were told this model is stationary when all of the roots of $\phi_p(B)$ are outside the unit circle, but I don't remember hearing anything about when these are invertible. If I use the parallels from the MA(q) model, I might say the AR(p) model is always invertible. Is that correct? Thanks


title: How to compute confidence interval for number to treat in a logistic regression?

text: I have built a multivariate regression model where $\mathbf x$ is the variable vector. Now I am asked to calculate the confidence interval for the number to treat (NNT) at $\mathbf x=\mathbf {x_0}$ which in turn requires to build the confidence interval for $p=Pr(1|\mathbf x=\mathbf {x_0})$. I'm using JMP which, to my knowledge, shows only point value and standard error of the parameter estimates but not the complete Fisher information. So I'm temped to use $\hat{p}=Pr(1|\mathbf x=\mathbf {x_0},\mathbf {\hat{b}})$ where $\mathbf {\hat{b}}$ is the vector of the ML parameter estimates, and then calculate $\hat{p}(1-\hat{p})/n$, but I doubt it is correct. Which is the simple way to do?


title: Kaplan-Meier statistic with no censoring reduces to empirical survival function

text: I just started taking survival analysis class and I'm stumped on this question. We need to show that when there are no censored observations $\hat{S}(t)=\prod_t_{(i)\le t}(n_i-d_i)/n_i$ equals to the empirical survival function $S_n(t)=\\# \\{ t_{(i)} \ge t \\}/n $ It looks pretty obvious to me so I'm not sure how to approach this.


title: Fit statistics and log-likelihood after multiple imputation estimation

text: I am working with four survey datasets, two of them without a certain variable which I impute using multiple imputation techniques. Afterwards, I use the complete data to estimate a set of logit models for which -to my understanding- there is no substantive interpretation for their log- likelihood. This is an important issue when you need to evaluate the relative goodness of fit of the models, since most fit statistics use the log- likelihood in their formulations: BIC, AIC, most pseudo R-squared, Deviance, etc... Is there a way to estimate a substantive log-likelihood for a logit model after multiple imputation? If not, which fit statistics would you recommend instead? PS: I'm working mainly in Stata, but also in R.


title: Computing probability of completing a task composed of independent events

text: This is a general question. I have a task that is composed of 3 independent events: A, B, and C. All are mutually exclusive and don't happen at the same time. So first A, then B, then C. I know the probability of completing each event with respect to time, I have the pdfs. How can I calculate the probability of completing the task at hand with respect to time?


title: I have a statistic, how do I calculate its distribution?

text: I am comparing images for correlation. The images are all correlated, but I would like to determine when one pair is much more highly correlated, relative to another pair. I am using as a statistic the difference in the logs of p-values produced by a Spearman rank correlation test. If this is difference is small, there is a normal degree of correlation. If this difference is big, there is a higher degree of correlation in one pair than the other pair. Let $A$, $B$, and $C$ be images. I would like to determine a distribution and calculate a p-value for $ x = \vert \log(\mbox{SpearmanRankCorr}(A,B)) - \log(\mbox{SpearmanRankCorr}(A,C)) \vert $


title: Danger of setting all initial weights to zero in Backpropagation

text: Why is it dangerous to initialize weights with zeros? Is there any simple example that demonstrates it?


title: Definition of "precision" of an estimator

text: What is the exact definition of _precision_ of a statistical estimator? Any references to textbooks are appreciated.


title: Regression with only categorical variables

text: Is it possible to conduct a regression if all dependent and independent variables are categorical variables?


title: Finding a prediction interval for an lmer model via mcmcsamp versus simulate

text: I've created a model with `lme4`'s `lmer` and wanted to create a prediction interval around my model fit. I figured I could do it Bayesian-style and simulate from the model to do this. My model looks like:               b3 <- lmer (dollars ~ 1 + I(refund/days) + (1 | month) + (1 | regime), data=elect5)      and my first attempt to graph this was with:               plot (elect5$dollars, type="o", ylim=c(200, 1000), ylab="Y", xlab="X")     for (i in 1:5000) lines (simulate (b3), col=rgb (0, 0, 0, 0.01))     lines (fitted (b3), col="green")     points (elect5$dollars, col="red")     abline (h=fixef (b3)[1], lty=3, col="blue")      which results in the graph (red points actual data, green line fitted model, gray "prediction interval"): ![Graph based on using simulate](http://i.stack.imgur.com/dBEZF.png) which doesn't look so good. It looks like it's simply a variance around the fixed effect intercept of the model. OK, so then I do things the "hard way" and code up a function to implement the model:               doB3mc <- function (x, f, r) f[1] + (f[2] * (x$refund / x$days)) + (r[x$month]) + (r[as.numeric (x$regime) + 12])      And then fill a matrix with simulations based on `mcmcsamp`, which samples the coefficients (not done the most efficient way, I'm sure):               a <- mcmcsamp (b3, 5000, saveb=T)     res <- matrix (NA, nrow=5000, ncol=49)     for (i in 1:5000) for (j in 1:49) res[i, j] <- doB3mc (elect5[j,], a@fixef[,i], a@ranef[,i])      And also graph:               plot (elect5$dollars, type="o", ylim=c(200, 1000), ylab="Y", xlab="X")     for (i in 1:5000) lines (res[i,], col=rgb (0, 0, 0, 0.01))     lines (fitted (b3), col="green")     points (elect5$dollars, col="red")     abline (h=fixef (b3)[1], lty=3, col="blue")      Resulting in the graph: ![Graph based on mcmcsamp](http://i.stack.imgur.com/y8V5k.png) Which looks more like what I'd expect, except it's biased: at values below the fixed effect intercept, the model is low in the prediction interval, and vice versa. So my questions are: 1) why does the first option, with `simulate` work strangely, 2) why does the second option look biased (fitted line not centered in the interval), and most importantly 3) is either of these the proper method for doing what I set out to do?


title: Comparing small and large samples from the same target population

text: Is it valid (reliable etc.) to compare a small sample with a large sample (which is made up of several small samples) when both samples are from the same population? There are 3 firms. Each firm specialises in the manufacture of a particular product, hence it is classified in that way. Firm A manufactures common cold medicines, Firm B manufactures arthritis medicines and Firm C manufactures heartburn medicines. The common element linking the three firms is that they are all medicine manufacturers. Now, Firm A and C only produce generic medicines and Firm B only produces branded medicines. I have a sample of 80 workers from each firm and I want to test how satisfied they are in regard to their work-life balance. I am using correlation analysis etc. I am comparing Firm A with B and C but I also want to compare Firm A and C with B. (Firm A and C are grouped together because they are generic medicine providers). Is comparing a sample of 80 peoples (Firm B) with a sample of 160 people (Firm A and C) any reason for statistical concern?


title: On testing for heteroskedascity, some questions on the White test & Breusch-Pagan test

text: 1. Isn't the test statistic for both tests identical? The only difference I see is that the alternative hypothesis is different. Since you don't need to know the function h in the Breusch-Pagan test, what exactly is the difference between the two tests?    2. When would you use the White test, and when would you use the Breusch-Pagan test instead? What's the most important criteria to consider here when deciding which test to use?


title: Ways to determine if experience or recent practice time is more significant in ranking?

text: **First a little background:** I'm a college paintball coach and I have two populations of players: I. Players that have a lot of experience/skill prior to joining the team e.g. playing for amateur non college teams but who don't often play together for various reasons. II. Players that are not very experienced/skilled but who have a lot of recent practice time playing together with other players in this population. I've also recently read this article by Michael Lewis about Shane Battier. In summary, it recounts how he has very low or non-existent traditional basketball "stats" but he has a very high "plus/minus" wherein whenever he is on the court, he limits the scoring of the other team. This leads me into the questions. **The questions:**   1. What is the best way to quantify (from a statistical and/or probabilistic sense) if experience/skill or recent playing experience is a better predictor of wins?   2. How do I identify players that have poor individual skills yet make large contributions to winning matches when playing with other players? **What I have so far:** a. Several individual skills tests results e.g. can you hit targets while running repeatedly b. A record of 1 vs 1 encounters that allows me to say things like: Player A has beat Player B 66% of the time in 1 v 1 matches c. A record of n vs n encounters (with n being up to 5 inclusive) Each weekly practice gives me additional data to use as well. My thinking is that after enough 1 v 1 matches each player would have played every other player multiple times. This would allow me to rank everyone at an individual level (e.g. using ELO). Next, I would take the ELO rankings and combine them and see if individual skills (data from group a and b) correlate to high number of wins (data from group c). I'm also not quite sure on the best way to combine the ELO rankings e.g. sum them, average them, take the min, take the max etc. Is there a way to structure the n v n matches in such a way that this becomes easier? e.g. if round robin every possible 3 v 3 combination will this increase the sample size and give me better eventual outcomes? As a side benefit, as mentioned above, if I could find players that have low individual skills scores and/or low ELO rankings but somehow have a high correlation with wins when playing on n v n matches (with n being greater than 1). If there is any more detail needed, please let me know. Thanks in advance. **UPDATE 2012-02-11** So we're now a couple weeks in of tracking how each player does in 1v1 and 3v3 match ups. I've created two ELO rankings based on each category (one for 1v1 and one for 3v3). I used the same formula and weight as this page. For the 3v3, I took the average in the same manner that the link does for doubles. For the 1v1 we have on average 6 match ups and for the 3v3 we have on average 9 so the ELO rankings are probably a little premature especially since we started everyone at 1500 for each category. Here is a scatter plot showing the two categories where each point is one player: ![scatterplot](http://i.imgur.com/5Ack2.png) My next question would be thoughts on how best to calculate the 3v3 ELO e.g. is it better to take the average, max, min, average of min and max etc. I'm keeping track of how accurate the ELO is for predicting the winner as well so hopefully the data will answer this question as well.


title: decision tree using user defined split function in rpart: No splits returned when tree is run

text: I've written a user defined splitting function to use with rpart, its returning a 'vector of goodness', but the tree that is returned never has any splits, just one node. using the anova method on the same data provides a tree with many nodes. I've been interrogating the code, but rpart feels like a black box so far, I can't figure out why there wouldn't be any splits or how to trouble shoot the problem. Not sure whether this belongs here or @stackoverflow. I'll include my function               distance <- function(y, wt, x, parms, continuous) {         # return a vector of goodness using a distance measure         # find distance between treatment and control in left and right,         # gain = difference between the two, normalized               if (missing(parms) || is.null(parms['treated']))             stop("requires vector indicating treatment")         treated <- which(parms['treated'] == max(parms['treated']))         #treated <- which(parms$treated == max(parms$treated))              goodness <- numeric(length(x) - 1)              for (i in 1:(length(x)-1)) {             left <- which(x >= x[i])             left.data <- y[left]             right.data <- y[-left]                  left.distance <- dist(rbind(left.data[treated], left.data[-treated]))              right.distance <- dist(rbind(right.data[treated], right.data[-treated]))                    treated.frac <- length(left.data[treated])/length(left)                  normalization <- gini(c(treated.frac, 1-treated.frac)) * left.distance             + treated.frac * gini(left.data[treated])             + (1 - treated.frac) * gini(left.data[-treated]) + .5                  goodness[i] <- (left.distance - right.distance)/normalization         }         list(goodness = goodness, direction=rep(1, length(x) - 1))     }


title: How to distinguish between periodic and random impulse?

text: I am working for two signals. One is a dataset with 10 equally spaced impulses. Another is a dataset with randomly spaced impulses of identical length. The higher order moments, such as kurtosis, skewness or standard deviation are not sensitive to the two signals. Who can help look for an indicator to identify it?


title: Warning:"NAs introduced by coercion" in MICE with unique ID

text: I am having a problem using MICE, where it generates the following warning:               Warning message:     In var(data[, j], na.rm = TRUE) : NAs introduced by coercion      This seems to be related to the presence of a unique alphanumeric ID variable in the data. For example:               require(mice)          imp <- mice(nhanes,maxit=0)     # no problem          nhanes$id <- 1:25     imp <- mice(nhanes,exclude="id",maxit=0)     # no problem          nhanes$id <- LETTERS[1:25]     imp <- mice(nhanes,exclude="id",maxit=0)      The latter generates the warning:               In var(data[, j], na.rm = TRUE) : NAs introduced by coercion      The predictor matrix is > imp$pred                   age bmi hyp chl id     age   0   0   0   0  0     bmi   1   0   1   1  0     hyp   1   1   0   1  0     chl   1   1   1   0  0     id    0   0   0   0  0      Note that here `id` is not used to impute anything. So it seems strange that a warning would be generated at all, in this simple case. Even this generates the same warning:               pred <- imp$pred          imp <- mice(       data=nhanes,       m = 5,         maxit=5,       imputationMethod = c(         "norm",  # age         "norm",  # bmi         "norm",  # hyp         "norm",  # chl         "")    ,  # id       predictorMatrix = pred     )      where id is not specified to be imputed or to be used for imputation in the predictor matrix, and the imputation method is "" So my questions are: >   1. Why is this warning generated ? >   2. How can I prevent this warning, without removing the ID variable from > the data. Of course I don't use the ID variable for imputing other variables > - it's just a bit inconvenient to remove it. >   3. If it can't be prevented, can I safely ignore the warning and use the > predictor matrix for imputing my data (after setting the id column to zeros) > ? >


title: Comparing divergence from uniform distributions with differing supports (discrete)

text: Imagine we have a potentially biased coin and a potentially biased six-sided die and we want to know which is more biased than the other. Firstly, is this a reasonable goal? Could it make sense to say, for example, that the coin is more biased than the die? Secondly, assuming the goal is reasonable, can we use the KL-divergence from the natural uniform (unbiased) distribution to the observed distribution as a measure of bias in such a setup? If not, what might we use instead? Lastly, if using KL-divergence is reasonable, can we obtain any meaning from the degree of difference between the divergences? For example, if the coin diverges by 2 and the die diverges by 1, can we say the coin is twice as biased as the die? The coin is to be compared with the uniform distribution over two outcomes while the die is to be compared with a uniform distribution over six outcomes. Clearly, if both discrete random variables were over the same number of possible outcomes, a direct comparison of the KL-divergences would be reasonable. But does this remain true when the number of the outcomes differs?


title: Back transformation of an MLR model

text: I've obtained a multiple linear regression model in the form $$ \mathrm{log}(Y) = \beta_0 + \beta_1x_1 + \dots + \beta_4x_4 + \beta_5x_1x_2 + \dots + \beta_{10}x_3x_4 + \beta_{11}x_1^2 + \dots + \beta_{14}x_4^2 $$ from the model I've obtained some contour plots that, for a more immediate interpretation, I've back transformed in original units (using $10^{\tilde{y}}$). Now my question is: do I need to back transform also the model parameters, the confidence intervals of the parameters themselves and the interval of prediction for new responses? How? I think that would be nice to back transform at least the RMSEP and RMSE in order to obtain an estimate in original units of the mean error associated with the model in prediction and in fitting. How? Thanks.


title: Who decided that p value should be less than 0.05 and why?

text: In many scientific fields, a p value less than 0.05 is considered significant, pointing to real difference in reality, whereas anything above it cannot be accepted (i.e. null hypothesis cannot be rejected, which means your results cannot be accepted for publication), From a perspective of pure statistics, this threshold of p value is completely arbitrary. Why have the scientific community settled on this value, and is willing to reject the null hypothesis only for for anything with a probability less than 5%?


title: Mann-Whitney U test with unequal sample sizes

text: I have two unequal groups (94 and 52) and want to run a Mann-Whitney U-test to see whether their scores on a measured variable differ. I see that it is OK to do with Kruskall-Wallis, does the same apply to Mann-Whitney?


title: Is use of reliability analysis always necessary in data handling?

text: My study involves studying the seasonal variations (premonsoon, monsoon, postmonsoon) in river at 5 different sampling points. The water quality study involved eight water parameters. To check the reliability of data, I had run reliability analysis and got Cronbach's value = 0.187..which is very very low. But my question is, is it necessary to run reliability analysis for every data? The data I am handling is completely scientific kind of data..which can show inconsistent results..yet is the standard reliability test necessary? If yes, please advice on what can be done about this as the shifting of data etc was of no use.


title: MC Standard Error for Metropolis Random-Walk Algorithm

text: I need to calculate some integrals using MCMC with Metropolis Random-Walk Algorithm. To decide the value I will accept as my integral, I will calculate 5 simulations of MCMC with size 2000, but I won't use the first 500 obtained values (Burn-in step). So, my doubt is: $$ SE = \sqrt{\frac{1}{N} \sum_{i=1}^N (g(X_i) - \bar{\mu}_N)^2}$$ where $\bar{\mu}_N$ is the estimative of $E(g(X))$ and $x_i \sim f$ for some distribution with pdf $f$ 1) Should I use the same generated sample $X=(x_1, ... , x_N)$ that I use to calculate $\bar{\mu}_N$ to calculate the standard error? 2) Should I consider the burn in sample to calculate the $SE$? That is, should I use $N=1500$ or $2000$? 3) I know that my integral is something between 0.98 and 0.97 and I obtained, in one simulation, the following: Integral:0.985333333333 $SE$:0.120214622886 Is that $SE$ correct? I mean, should it be smaller? Thanks in advance!


title: Closed form expression for the distribution of the sample kurtosis of Gaussian distribution

text: Is there a closed-form expression for the distribution of the Sample Kurtosis of data sampled from Gaussian distribution? i.e., $P(\hat{K}<a)$ where $\hat{K}$ is the sample kurtosis.


title: Accounting for Violation of Parallel Trend Assumption in Diff-in-Diff with Propensity-Score Matching

text: Objective: I want to test the effect of a regulatory change using a classical pre-post/treatment-control DiD design (Y = POST + TREAT + POSTxTREAT + e). Problem: Treatment & control obs. are from different countries. Although these countries are similar with respect to many relevant aspects, a visual check of the parallel trend assumption for the outcome variable in the pre-period (only 3 years) indicates different time trends. Sample details: time period is 2002-2007, event happened in 2005. Treatment obs. 200/year (1200 total), control obs. ~ 100 obs./year (600 total). Obs. over time belong to the same subjects (300 subjects). Questions:   1. Does it make sense to use propensity-score matching based on determinants of the outcome variable in the pre-period to create a more homogenous treatment/control sample (and hopefully get rid of the different trends)?    2. What would be the advantage over simply controlling for the determinants used in the PSM in the DiD regression? (e.g., Y = POST + TREAT + POSTxTREAT + CONTROLS + e)   3. I thought of using one-on-one matching with a caliper. However, I have more treatment obs. than control obs. – are there any better ways to avoid losing too many of my obs.? Thanks! PS: I'm using Stata.


title: Questions on likelihood analysis

text: Whilst studying likelihood methodologies, I've come across some results that I haven't been able to work out.   1. If $X$ and $Y$ are Poisson with means $\mu_{X}$ and $\mu_{Y}$, then the conditional distribution of $X$ given $X+Y$ is binomial with parameters $n = X+Y$ and $\pi = \frac{\mu_{X}}{mu_{X}+mu_{Y}}$.   2. Assuming $y_{1}, \ldots, y_{n}$ are independent exponential outcomes with mean $\mu_{1}, \ldots, \mu_{n}$. Given $\text{log}(\mu_{i}) = \beta_{0}+\beta_{1}x_{i}$ and $\sum x_{i} = 0$, the profile likelihood of $\beta_{1}$ is $$ \text{log}(L_{p}(\beta_{1})) = -n \log \left(\sum_{i}y_{i}e^{-\beta_{1}x_{i}}\right)$$ I would include my attempts at showing these results, but I'm completely lost here. Perhaps somebody could provide some assistance.


title: R - Post hoc comparisons in 2x2 ANOVA unequal sample sizes

text: I have a 2x2 between-subject design with (slightly) unequal cell sizes. Hence I chose to use the `anova(lm(...))` function so that Type 1 SS are used: `anova(lm(prob ~ utility * imagination, data = e2data))` `prob`is the numerical DV and `utility` and `imagination` are factors with 2 levels each. The ANOVA yields **no significant main or interaction effects**. Then, I went on to visualise the cell means including errorbars:               library(reshape)     mystats <- function(x) round(c(M = mean(x), SD = sd(x), n = length(x), SE = sd(x)/sqrt(length(x))), 2)     prob.means1 <- data.frame(cast(e2data, utility + imagination ~ ., value="prob", mystats))      `mystats` calculates the mean, standard deviation, number of observations and standard error of the mean. This results in the following data.frame, which is the basis for the plot below:               > prob.means1       utility imagination    M   SD   n   SE     1    gain      absent 7.18 1.23 100 0.12     2    gain     present 7.49 1.34  95 0.14     3    loss      absent 7.10 1.47 103 0.14     4    loss     present 7.19 1.49  95 0.15      The plot:               (plot1 <- ggplot(prob.means1, aes(x = utility, y = M,                                       group = imagination,                                       color = imagination)) +        geom_point(size = 3.5) +        geom_line() +        geom_errorbar(aes(ymax=M+SE, ymin=M-SE, width=0.05)) +       scale_y_continuous(name="Probability estimate", limits=c(6,8))     )      Note that I set the errorbars to range from M-SE to M+SE. ![2x2 plot](http://i.stack.imgur.com/sSybJ.png) **In the gain condition (left part of the plot) the error bars do not overlap. Does this allow the conclusion that the difference between these two cells is significant? How do I run this post-hoc comparison appropriately (between- subject, unequal sample sizes)?** I tried `Tukey.HSD()` but it doesn't work on a `lm` object. Also, I tried to use `pairwise.t.test` but I couldn't figure out how to use it with more than 1 factor.


title: Similarity between different length vectors containing related items

text: I have a vector (V1) with which I need to calculate the similarity of other vectors (ex V2,V3 ... ) which may be of different lengths. The different angle here is that the elements inside the vectors are itself similar to each other.               V1 = c("a","b","c")     V2 = c("a","d","e","f","g")     V3 = c("b","c","f")      The elements are similar to each other for example:                     a      b        c       d       e       f      g     a   1.00    0.18    0.01    0.96    0.12    0.46    0.73     b           1.00    0.07    0.36    0.13    0.47    0.92     c                   1.00    0.88    0.62    0.65    0.31     d                           1.00    0.86    0.96    0.55     e                                   1.00    0.25    0.91     f                                           1.00    0.13     g                                                   1.00      The usual similarity methods like cosine similarity, correlation do not make use of the similarity between the elements of the vectors. Any suggestions to how should I calculate the vector similarity?? I am using R.


title: T-Testing, test if one average exceeds another at the 5% level

text: I struggle with following question: > It is now considered that house prices are rising again following the recent > slump. Given that in 2008 the average price of a standard 3 bedroom house > was €215,000. A sample of 25 similar houses in 2013 gave an average price of > €275,000 with a standard deviation of €50,000. Test to determine if the > average price for the houses in 2013 exceeds the average price in 2008 at > the 5% level. Yes\No answer only, will receive zero marks. NOTE: NO dataset is provided. This is an assignment question which of course I don't want you to do it all for me!! I'm bit confused as can't figure out what equation/formula to use. As reference I'm using my slides notes and Statistics boook (James McClave, 2014) 12th edition. I'm using this formula: http://www.nursingplanet.com/biostatistics/inferential_statistics_t-test.html Student t-testing calculating t-statistic for comparing means. Any help hints/tips would be appreciated! Thanks!


title: Standard Deviation of Random effect is 0?

text: I have a model with two random effects               > lmer(TotalPayoff~Type+Game+PgvnD*Asym+(1|Subject)+(1|Pairing),REML=FALSE,data=table)-         >m1      each pairing contains 2 subjects and each Subject is in 2 different pairings with different subjects. For some reason the subject random effect has SD =0 I do not understand this? Is it something wrong with how I structured my random effects?               > m1     Linear mixed model fit by maximum likelihood      Formula: TotalPayoff ~ Type + Game + PgvnD * Asym + (1 | Pairing) + (1 |      Subject)         Data: table        AIC  BIC logLik deviance REMLdev      1014 1038 -497.8    995.6   964.4     Random effects:      Groups   Name        Variance Std.Dev.      Subject  (Intercept)   0.000   0.0000       Pairing  (Intercept) 716.102  26.7601       Residual              89.364   9.4532      Number of obs: 113, groups: Subject, 73; Pairing, 61          Fixed effects:                 Estimate Std. Error t value     (Intercept)   58.428     12.255   4.768     Type           7.926      2.852   2.779     Game          15.374      7.147   2.151     PgvnD         -8.466      7.554  -1.121     Asym         -12.167      7.583  -1.604     PgvnD:Asym    26.618      9.710   2.741          Correlation of Fixed Effects:                (Intr) Type   Game   PgvnD  Asym       Type       -0.335                                 Game       -0.835  0.009                          PgvnD      -0.098 -0.038 -0.010                   Asym       -0.330  0.081 -0.015  0.189            PgvnD:Asym  0.189 -0.267 -0.011 -0.766 -0.328


title: Examples for commercial applications of data mining in pharmacy

text: I've always heard that data mining and machine learning tools/techniques are heavily used in the pharmacy sector and biology but I have never heard of companies that offer commercial applications of data ming for pharmacy and biology. Once I heard my professor saying that if you really want to know the state-of-the-art in data mining then look at the pharmacy sector. But again I can't find any companies that apply them. My questions: 1- Can you please tell me some commercial applications and usage of data mining in the pharmacy and biology sector? 2- What companies that offer such applications for the pharmacy and biology sector?


title: How to calculate the regression variance for a GLS model?

text: I need to calculate the regression variance ($\sigma^2$) in order to estimate both the confidence intervals and the prediction intervals in a gls regression analysis. For the analysis, the covariance matrix ($V$) of the response variable ($y$) is known in advance, and so I use it directly as the weighting matrix (=$V^{-1}$) in the gls regression analysis. The regression variance is a weighted sum of the residual error: $\sigma^2 = \frac{ (Y – X\beta)^T C^{-1} (Y – X\beta)}{n – p}$ My question/problem is how to determine the weighting matrix $C^{-1}$? $C$ cannot be set equal to $V$ since (according to the above equation) $C$ must be dimensionless while $V$ has the same units as $\sigma^2$. Based on my reading of the literature and available texts, it seems that $C$ is the correlation matrix and is a scaled or normalized form of the covariance matrix $V$. i.e., $V = Var(\epsilon^2) = \sigma^2 C$. But my problem is that $\sigma^2$ is not yet known, and so I need another way find $C$ from $V$. R functions such as gls() will compute the regression variance (if I knew how gls() does this, it would answer my question). However I cannot use gls() in this case since I am specifying a user-defined covariance (weighting) matrix, and gls() only accepts a limited set of specific correlation structures. In fact a possible solution can be found in this earlier post where an equation for the SEE (or sigma2) for a GLS regression was cited : GLS calc of SEE: sqrt( sum( ( residuals from linear model) ^ 2 * glsWeight ) ) / sum( glsWeight ) * length( glsWeight ) / residualDegreeFreedom ) However I am unable to ascertain the validity of this equation and cannot find its source reference.


title: Quantifying Potential for Violence/Espionage

text: I am planning a seven-question Likert-type survey with a "1 to 10" scale. The questions I am asking in the survey are to assess an individual's risk factors (how much risk they present to an organization). The risk factors are compounding, meaning the more risk factors an individual exhibits, the greater risk the individual represents; so the higher they rate an answer on the response scale, the greater the indication for that questions. Rather than comparing the Likert-items to each other as is done in a standard Likert survey, the questions are normalized per respondent and then measured against the population. I plan to take the geometric mean of each respondent's answers (because, while the data is ordinal with an indefinite distance, each respondent infers their own distance between ordinals based on their background, and a median does not account for compounding or outliers) and then create a distribution of all respondents. What I want to do is create a Pareto cut-off so that the lowest 20% of scores in the range (not the population) will be excluded. To do so, I am adding 20% of the range to the lowest response in the population:               ((xMax - xMin) * .2) + xMin = cut-off      What I am finding is that when calculating the right-tail P-Value of the normal cumulative distribution (using computer generated random test data) I am always right around .8, typically somewhere between .79 and .82 for P-Value.   * Is this just an anomaly due to using a weak pseudo-random number generator?   * If I take the lowest score in the population and add 20% of the range, it is not the same as finding the 20th percentile of the population, right?   * Is there a better way to do this? Thanks! Justin


title: Generalized linear model: link function Power(-1)

text: During study our of statistics in my psychology coursework, we had to teach ourselves how to use generalized linear models in SPSS (only basic knowledge). For an exam we may also use generalized linear models and I want to try it, but I have two problems:   1. According to the distribution of the dependent variable, a normal distribution would fit our data best. A fellow student said that this is not logical, because our dependent variable is "time required to complete a task", so it cannot have negative values. Our lecturer wrote on his slides that a variable needs to be able to range from $-∞$ to $+∞$ to use the normal distribution.  In another post in these sides, I read that the normal distribution creates a predicted variable which fits these characteristics, so the original dependent variable doesn't need to fulfill those characteristics. Is that true?   2. As I can't use the normal distribution, a Gamma distribution would fit best. According to SPSS's documentation, the canonical link function is Power(-1). I don't quite understand what I have to put in the box "Power" in SPSS. Do I have to fill in 0,05 if I want to work with a Power of 0,95 (so Power(-1) would mean 1-0,95) or do I have to fill in -1 as value of power (which seems quite unusual to me).


title: How do I test whether two regressions lines cross at y = 0?

text: I have measurements for daily space heating vs daily mean outdoor temperature for two different control strategies. The data is shown here: ![enter image description here](http://i.stack.imgur.com/c4L6Y.png) I have also performed linear regression, of the form:               lm(Energy ~ Control * MeanOutdoorTemp)      This yields four coefficients:               Coefficients:                                Estimate Std. Error t value Pr(>|t|)         (Intercept)                170.8293     4.2083  40.594  < 2e-16 ***     ControlMPC                 -30.7044     4.9025  -6.263 1.38e-07 ***     MeanOutdoorTemp             -6.2924     1.6466  -3.821 0.000413 ***     ControlMPC:MeanOutdoorTemp   0.8211     1.7162   0.478 0.634709          I'd like to test whether these two regression lines cross at zero energy. What would be the statistically correct way to do that?


title: Can I use logistic regression when my IV is dichotomous but my DV is continuous?

text: I have a longitudinal data set in which my predictor/Independent Variable is dichotomous and was assessed at age 4 years, and my outcome/Dependent Variable (depressive symptoms) was assessed at 16 years. I would rather not arbitrarily dichotomize the DV, and I would prefer to avoid non-parametric stats so that I could include some covariates. Naively, I would think that one could use a dichotomous IV but continuous DV in logistic regression, though I would assume that technically I would be using depressive symptoms to predict my IV? Is there any precedence for this in the literature (or alternative analyses to consider)? Thanks.


title: MLE2 (in R): I can't find the right starting values

text: I am trying to fit several MLEs on a small dataset. I have managed to get it done for all my potential models except for this last one. This last one is actually the model that I expected to be the best model, because it was the best model when I did the whole thing using GLS. To start, here is the data:               dat$rate <- c(8,1,17,34,8,30,8,15,17,4,29,12,29,12,12,18,24,20,10,2,4)     dat$pred <- c(2,1,2,2,1,3,1,3,2,1,1,2,3,2,4,4,4,4,3,1,2)     dat$prey <- c(10,10,20,60,20,40,30,20,30,40,50,40,60,50,20,30,40,50,20,10,10)     dat$mass <- c(6.885000,5.942000,7.199000,11.285500,7.263000,6.571667,  4.583000,4.949000,6.220000,2.730000,13.330000,6.074000,7.053667,6.626000, 3.291500,6.711750,4.530750,7.076750,6.644333,2.861000,3.793500)      Here is the function that contains the equation:               CM = function(N0,a,h,tt,P,m, S, al) {     a = a/(1+m*(P-1)+a*h*m*(P-1)*N0)     h=h*S^al     N0 - lambertW(a*h*N0*exp(-a*(P*tt-h*N0)))/(a*h)     }      And here the model:               model = function(eaten, initial, attack, handling, time.exp, pred.no, interf, mass, allometric) {     if (attack < 0 || handling < 0) return(NA)     prop.exp = CM(N0=initial, a = attack, h = handling, tt = time.exp, P = pred.no, m=interf, S= mass, al= allometric)/initial     likelihoods = dbinom(eaten, prob = prop.exp, size = initial, log = TRUE)     negLL = -sum(likelihoods)     return(negLL)     }      And finally the code to fit the model               cm3 = mle2(model,            start = list(attack = .02, handling = .02, interf=0.02, allometric = -0.2),            data = list(initial = dat$prey,                            eaten = dat$rate,                        pred.no = dat$pred,                            mass = dat$mass            ),            fixed = list(time.exp = 1)     )      I have tried many different starting values and they all ended up with error messages like:                non-finite finite-difference value [2]      or                initial value in 'vmmin' is not finite      I guess that either "attack" or "handling" returns an invalid value but I am not sure how to fix this. Any ideas or suggestion are really appreciated. Thanks in advanced, Robbie


title: Is this a scale or ordinal variable

text: I have a questionnaire distributed to the users with the following question from which users can select an answer Q: What was the project Cost Performance Index               Option 1 : 0 < CPI < .5          Option 2 : .5 < CPI < .75          Option 3 : .75 < CPI < 1          Option 4 : 1 < CPI < 1.5          Option 5 : 1.5 < CPI       Each of the options have a weight assigned to them starting from 1 to 5. Is this variable continuous or ordinal?


title: Gamma linear regression, errors with two factor interaction. None with three factor

text: First off, I am aware I have not enough data to be performing a regression with this number of parameters, lets overlook that for now - I'm interested only in why I get an error with one model, and not with the other: I try to do a regression as follows:               glm("y~(x1+x2+x3+...x8)^2", family=Gamma(link=identity), data=data)      and I get:               Error: no valid set of coefficients has been found: please supply starting values      I try to do:               glm("y~(x1+x2+x3+...x8)^3", family=Gamma(link=identity), data=data)      and I get:               Warning message:     glm.fit: algorithm did not converge      but no errors, as I understand it the first error is complaining that I dont have enough data to fit the model by itself, so starting coefficients are required. Why do I get this error when trying to utilise less predictors, and not when using more?


title: Duality between acceptance region of testing and confidence region in terms of optimality?

text: The duality theorem between acceptance region of testing and confidence region is about their validity (i.e. satisfying the significance level of the former and the confidence level of the latter sum up to be 1). For example, this is what is given in Casella's Statistical Inference. I was wondering if there is some duality result about their optimality? By "optimality":   * for confidence interval, it can mean the length of the interval is minimum, or something else;    * for acceptance region of testing, it can mean the power of the testing has uniformly maximum power, or the length of the acceptance region is minimum. Thanks and regards!


title: Contraindication for STL decompostion

text: Trying to "understand" a time series' patterns it is intuitively tempting to use STL decomposition as the concept of distinguishing between trend, season and the rest makes sense. But my experience tells me that no static algorithm will lead under all circumstances to useful results. So my general question/s is/are when should you not apply STL decomposition and if you do, what observation in the STL result might in you experience indicate a faulty/useless decomposition? Like you wouldn't blindly trust a correlation analysis of two variables without having a look at a scatter plot, b/c outliers might lead to a high correlation coefficient indicating a non-existant relation. I'm a newbie in this area, so a more extensive answer would be great.


title: Interpretation of mixed model coefficients

text: I ran a mixed model regression and obtained the following output. Because I log-transformed the concentration variable, I interpreted the coefficient as follows: the mean concentration is reduced $1.5$ percent each year. The years run from 1993-2010. Can say also state that the mean concentration is reduced by $26\%$ ($0.01509\cdot 17$) over the $17$ year period? Why or why not? Why do most people generally state $1$ unit increase in $X$ corresponds to a certain number change in $Y$ rather than infer change over the entire period?               baggerTrend <-  lme(Log.Qconc   ~   Yearf,  random=(~Yearf|MineID), Bagger)           Random  effects:                             Formula:    ~Yearf  |   MineID               Structure:  General positive-definite,  Log-Cholesky    parametrization                    StdDev    Corr                     (Intercept) 0.209172851 (Intr)                   Yearf       0.000374953 -0.478                   Residual    0.785367538                                Fixed effects:  Log.Qconc ~ Year                                        Value   Std. Error    DF      t-value       pvalue     Intercept   33.71122    6.3901     2762      5.275538       0     Year        -0.01509    0.003193   2762     -4.724678       0


title: Is a chi-sq test for independence appropriate for categorical time series?

text: I'm analyzing a very large, but not finite, set of survey data. I've assembled a time series crossing gender (2 categories, if there was uncertainty) with age (in 6 categories) and wish to test whether the demographic representation has changed over time. Using the chi-squared test for independence, I understand that failing to reject the null hypothesis means the gender-age distribution between two data points is unrelated statistically, and I can confirm this at points where we know a substantial change was made (for example, an increase in total surveys taken). There is zero expected overlap of respondents in two different time periods (a control in the survey process itself). My question is, since this is a repeated test (23 adjacent month-to-month pairs of distributions from 24 months of observed distributions), am I amplifying the type II error, and if so is there a better test?


title: Expectation of weakly stationary process

text: Hi all I have a question. Suppose $$X(t) = A \cos(t) + B \sin(t)$$ is a weakly stationary process, where $A$ and $B$ are random variables. Why is it that the mean of $A$ and $B$ are necessarily 0? As I understand it, weak stationarity only requires that the $E(X) = \mathrm{constant}$.


title: Estimating a survival probability in R

text: Based on a sample of $n$ survival times, I would like to estimate the probability of surviving time $t$, for some specific $t$, using the Kaplan- Meier estimator. Is it possible to do this in `R`? Please, note that $t$ is not necessarily an event time.


title: Kendall's tau for holdouts low and not significant - Conjoint

text: I have done my Conjoint Analysis (fractional factorial design) but when it comes to validating the model, it shows a Kendall's tau for Holdouts of 0.33 and not significant. But Pearson's R and Kendall's tau are high and significant which is supposed to show that the fit of my model is good... How do I interpret this result? When I remove 2 of my 4 holdouts, then it shows perfect Kendall's tau & significant...But I can't just take 2 holdouts out to 'make it work'... Thanks for your help. PS: Also, for my understanding, what regression or analysis is that conjoint command in SPSS running exactly?


title: What is the confidence interval calculated in a spectral density periodogram in R?

text: This question is similar to the one posed here: Testing significance of peaks in spectral density In that post, Pantera asks how to test whether a peak in a periodogram has a spike that is significantly different from the spikes generated by noise. My question is this, when using the spec.pgram() function in R, a confidence interval is calculated and plotted as a blue line -- what is the correct interpretation of this interval, and how can we extract what its numerical values are? I've seen it suggested that to find an answer to this question, we should inspect the code for the plot.spec() function and then observe the spec.ci() function code, and figure things out from there.               spec.ci <- function(spec.obj, coverage = 0.95) {           if (coverage < 0 || coverage >= 1)                stop("coverage probability out of range [0,1)")           tail <- (1 - coverage)           df <- spec.obj$df           upper.quantile <- 1 - tail * pchisq(df, df, lower.tail = FALSE)           lower.quantile <- tail * pchisq(df, df)           1/(qchisq(c(upper.quantile, lower.quantile), df)/df)       }      ![Plot generated by specpgram](http://i.stack.imgur.com/mN8x5.png) I suspect that the height of the CI is some indication of the amplitude for which 95% of the data are less than. But what is the "baseline" for that? Would it be possible to calculate and add horizontal lines to the plot indicating which peaks are significant? Any help would be greatly appreciated. For context, I am a biologist and these data are from imaging of rhythmic behavior in insects.


title: "Investigator intention" and thresholds/p-values

text: I'm reading John Kruschke's "Doing Bayesian Data Analysis" slides, but actually have a question about his interpretation of t-tests and/or the whole null-hypothesis significance testing framework. He argues that p-values are ill-defined because they depend on the investigator's intentions. In particular, he gives an example (pages 3-6) of two labs that collect identical data sets comparing two treatments. One lab commits to collect data from 12 subjects (6 per condition), while the other collects data for a fixed duration, which also happens to yield 12 subjects. According to the slides, the critical $t$-value for $p<0.05$ differs between these two data collection schemes: $t_{\textrm{crit}}=2.33$ for the former, but $t_{\textrm{crit}}=2.45$ for the latter! A blog post--which I now cannot find--suggested that the fixed-duration scenario has more degrees of freedom since they could have collected data from 11, 13, or any other number of subjects, while the fixed-N scenario, by definition, has $N=12$. Could someone please explain to me:   * Why the critical value would differ between these conditions?    * (Assuming it's an issue) How one would go about correcting/comparing for the effects of different stopping criteria? I know that setting the stopping criteria based on significance (e.g., sample until $p<0.05$) can inflate the chances of a Type I error, but that doesn't seem to be going on here, since neither stopping rule depends on the outcome of the analysis.


title: marginal effects (and confidence interval) for interaction variables (Stata)

text: I'm trying to calculate marginal effects (and their confidence intervals) for an interaction variable. I'm using Stata and have panel data (pooled cross-sectional time series). I am using the `xtpcse` command. For some reason the results calculated using Stata's `margins` command produce exactly the same dy/dx outcome for all of the different values at which the marginal effects are set. I'm clearly doing something wrong, but I don't know what it is. Any help appreciated. **Edit:**               xtpcse outcome  var1 var2 var3[lagged outcome] var4 var5 var6 var7 var8[time invariant] var9 var10 time var4*var2 [interaction] var1*var2 [interaction]          margins , dydx ( var1) at (var2=(-10 -5 0 5 10))          Output:      Average marginal effects                          Number of obs   =        468     Model VCE    : Panel-corrected          Expression   : Fitted values, predict()     dy/dx w.r.t. : var1          1._at        : var2     =      -10          2._at        : var2     =      -5          3._at        : var2     =       0          4._at        : var2     =       5          5._at        : var2     =       10          ------------------------------------------------------------------------------                  |            Delta-method                  |      dy/dx   Std. Err.      z    P>|z|     [95% Conf. Interval]     -------------+----------------------------------------------------------------     var1  |              _at |               1  |   .0058414   .0019627     2.98   0.003     .0019946    .0096882               2  |   .0058414   .0019627     2.98   0.003     .0019946    .0096882               3  |   .0058414   .0019627     2.98   0.003     .0019946    .0096882               4  |   .0058414   .0019627     2.98   0.003     .0019946    .0096882               5  |   .0058414   .0019627     2.98   0.003     .0019946    .0096882


title: Mixed-effect modeling with paired observations & bounded response variables

text: I am quite new in the field of mixed-effect modeling. For a beginner like me, I guess I combine several levels of complication in my analysis: paired observations & bounded response variables. I hope to be clear in my description and to get some useful comments, at least on some of the challenges I face. Thanks in advance. In a nutshell, my objective is to compare performance of several species distribution model types based on several performance indexes, across a set of 89 species. **Responses variables:** AUC metric (bounded between 0 and 1) or Boyce index (bounded between -1 and 1) **Fixed effect:** model_type (climatic_model, land_use_model) **Random effect:** 89 species, 10 replicates for each of the 89 species I guess one could do: `fit1<-lme(model_performance_value~model_type, random=~1|species, data=dat)` Using summary(fit1), I get:               Fixed effects: model_performance_value ~ model_type                              Value   Std.Error   DF  t-value p-value     (Intercept)             0.6362836 0.028055699 1690 22.67930       0     techniqueunbiased_boyce 0.1680262 0.007156959 1690 23.47732       0      (1) Is it fine to get the 0 p-values (became <0.001 using aov(fit1))? (2) Do you think the DFs are ok? (3) If I use the AUC as response variable it is limited to [0,1] and if I use the Boyce index, it is restricted to [-1,1]. Should I use other models with such response variables? **CONSIDERING PAIRED COMPARISONS:** The response variable of species A replicate 1 climatic_model is paired with the response variable of species A replicate 1 land_use_model. The replicates are paired between model types but within a given species. I made another variable called paired_id. This variable was written as "species"_"replicate", and takes 10 replicates x 89 species=890 values.               model_type species replicate paired_id     clim          A        1         A_1     land          A        1         A_1     clim          A        2         A_2     land          A        2         A_2     clim          B        1         B_1      Following this link (example in a much easier context) http://www2.compute.dtu.dk/~hmad/GLM/Slides_2012/week09/lect09.pdf, I did:               fit1<-lme(model_performance_value~model_type, random=~1|species/paired_id, data=subset)      I obtained a 0 p-value, and a DF=890. I am not sure it works this way, at all, though. I guess it is more complicated than this. Has anyone some tips for me? Any help much appreciated!


title: Fisher Information as measurement of curvature and hence information

text: (We may just assume single parameter model and single variable data for this discussion). I understand that the "observed information" evaluated at the MLE measures the curvature and hence gives a gauge of how much information our data is giving us. "Observed information" here refers to having already obtained a sample and then computing $-d^2l/d\theta^2$ ($l$ here is log-likelihood). I also understand the definition of Fisher Information, and I understand that $l$ is random before the data is drawn. However, in terms of an intuitive explanation, Wikipedia and various other sources say that Fisher Information gives a measure of curvature or still tells us information in the same way as I described for observed information. This, I disagree. When I evaluate $I(\theta)$ by average the curvature over all possible $l$, we must keep in mind that $\theta$ need not be the peak of $l$ - for some $l$, it is, but for other $l$, we may be way, way away from the peak, and the curvature that these scenerios contribute to the Fisher Information may lead to funny effects. Why would we care about the curvature at $\theta$ for these weird $l$ where $\theta$ is clearly not a good estimation? Basically, I am not very convinced on the intuitive explanation given for Fisher Information. Is the intuitive explanation commonly given a correct one?


title: How to rank predictor variables in logistic regression in R

text: I have created a logistic regression model in R. The dependent variable is binary whether or not a landscape transitioned to another cover class denoted by 1s and 0s, and the independent or predictor variables are some characteristics of the landscapes such as elevation, slope, distance to urban areas, etc. I need to illustrate the difference between predictor variables in terms of their importance or their effect on the dependent variable (land cover change). I was wondering if there is a way to rank and plot variable importance in logistic regression in R. Any help is greatly appreciated! Thanks, Azad


title: Converting ARMA models to infinite AR process in R

text: I'm trying to recreate some test statistics that use the infinite AR representation of normal ARMA models. I found out about the function ARMAtoMA but have not been able to find the same functionality for AR representations as well. Do you guys know of a package that would make extracting them easy? I was hopeful when I found the ARMAtoMA package, because it allows for such easy implementation in iterative loops, so I kind of dread having to implement this by hand. Found the incredibly easy solution in this forum post: http://www.talkstats.com/showthread.php/8578-A-question-about-ARMAtoMA You just have to switch up the PHI and THETA in the ARMAtoMA module to get the PIs of the infinite AR process. I guess I just didn't think this through correctly.


title: Expected number of duplicates (triplicates etc) when drawing with replacement

text: I have the following problem: I have 100 unique items (n), and I'm selecting 43 (m) of them one at a time (with replacement). I need to solve for the expected number of uniques (only selected once, k = 1), doubles (selected exactly twice k = 2), tripples (exactly k =3), quads etc... I've been able to find plenty of results on the probability of there being at least one double (birthday paradox), but not on the expected number of pairs in the population.


title: AIC, BIC, DIC, model selection criteria

text: I am trying to understand the difference between these parameters, and their application. Was hoping to get some correction/clarification to my statements. I have a training set and cross-validation set. I am using bayesian statistics, I realize AIC and BIC can't be used in this context. AIC and BIC assess how well the model fits the training set. The DIC assess how well the model will fit some future data (ie the cross-validation set). So if using MLE methods, should one compute AIC and BIC on the cross- validation set? Conversely, when using bayesian methods, is there no need to compute the DIC for the cross-validation set because it by definition is related to future data? Thanks Mark


title: Finding the variance when you know the population mean

text: So my basic understanding is that I may end up in one of the following situations.   * I have a sample mean and therefore calculate a sample variance   * I have a population mean sand calculate a population variance. Though in one book I have been reading it looks like teh numerator of the variance formula is using the population mean in the variance formula but subtracting it from a sample value. eg variance = (xi - population mean)^2 I don't get why if you knew the population mean why wouldn't you calculate the variance off the data you used to get the mean ? Why would you calculate the variance just based off a sample ? I realize there ,must be some explanation but with my limited understanding of this I am unsure of what it is . thanks