title: Statistical terminology: Is my wording correct?

text: I have written: > A predetermined level of significance ($\alpha$ – level) was set to assess > the null hypothesis based on a probability value ($p$ value). If the > $p$-value is greater than the $\alpha$-level, then the null hypothesis can > be accepted. If the $p$-value is less than or equal to the $\alpha$-level > then the null hypothesis is rejected. For this analysis the $\alpha$-level > was set at $0.025$ giving a $95%$ confidence of the result. > > For the $\alpha$-level of $0.025$, it is $95%$ certain that the true > difference between population mean values lies within the confidence > interval. The estimate for difference is the difference between the > population mean values as calculated from the sample data. > > The $p$ value was calculated as shown: > > $$\displaystyle p = \frac{\bar{x_1} - > \bar{x_2}}{\sqrt{\frac{\sigma_1^2}{n_1} + \frac{\sigma_2^2}{n_2}}} $$ > > using $n-2$ degrees of freedom. Is the terminology in that correct? Have I used "$p$ values" correctly there or should it be $t$ value? I'm pretty sure my formula is correct as I can't justify equal variances for both samples. Also, is my conclusion to my hypothesis test correct? Also, am I right in comparing my $p$ values to my $\alpha$ values or should that say something else?


title: $N(\theta,\theta)$: MLE for a Normal where mean=variance

text: $\newcommand{\nd}{\frac{n}{2}}$For an $n$-sample following a Normal$(\mu=\theta,\sigma^2=\theta)$, how do we find the mle? I can find the root of the score function $$ \theta=\frac{1\pm\sqrt{1-4\frac{s}{n}}}{2},s=\sum x_i^2, $$ but I don't see which one is the maximum.   I tried to substitute in the second derivative of the log-likelihood, without success. For the likelihood, with $x=(x_1,x_2,\ldots,x_n)$, $$ f(x) = (2\pi)^{-n/2} \theta^{-n/2} \exp\left( -\frac{1}{2\theta}\sum(x_i-\theta)^2\right), $$ then, with $s=\sum x_i^2$ and $t=\sum x_i$, $$ \ln f(x) = -\nd \ln(2\pi) -\nd\ln\theta-\frac{s}{2\theta}-t+\nd\theta, $$ so that $$ \partial_\theta \ln f(x) = -\nd\frac{1}{\theta}+\frac{s}{2\theta^2}+\nd, $$ and the roots are given by $$ \theta^2-\theta+\frac{s}{n}=0. $$ Also, $$ \partial_{\theta,\theta} \ln f(x) = \nd \frac{1}{\theta^2} - \frac{s}{\theta^3}. $$


title: determine the growth rate given a time series data

text: I have a data sets like this:               structure(list(Month = structure(c(14975, 15006, 15034, 15065,      15095, 15126, 15156, 15187, 15218, 15248, 15279, 15309, 15340,      15371, 15400, 15431, 15461, 15492, 15522, 15553, 15584, 15614,      15645, 15675, 15706, 15737, 15765), class = "Date"), TotalLogins = c(284697404L,      268944957L, 297847827L, 287150001L, 277779620L, 262275285L, 284271058L,      294965702L, 285132804L, 238847338L, 287527433L, 314483537L, 324823553L,      322896485L, 329044914L, 318228530L, 324395065L, 324988644L, 335464023L,      336269471L, 324063033L, 349017727L, 347193478L, 355561387L, 373885187L,      356774443L, 386372600L)), .Names = c("Month", "TotalLogins"), class = "data.frame", row.names = c(11L,      8L, 18L, 1L, 21L, 16L, 14L, 4L, 27L, 25L, 23L, 6L, 12L, 9L, 19L,      2L, 22L, 17L, 15L, 5L, 28L, 26L, 24L, 7L, 13L, 10L, 20L))          summary(lm(TotalLogins ~ Month, data=xx))          Call:     lm(formula = TotalLogins ~ Month, data = xx)          Residuals:           Min        1Q    Median        3Q       Max      -59901492  -4717787   1266082   9432626  26785529           Coefficients:                   Estimate Std. Error t value Pr(>|t|)         (Intercept) -1.674e+09  2.121e+08  -7.892 3.01e-08 ***     Month        1.294e+05  1.380e+04   9.377 1.16e-09 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Residual standard error: 17020000 on 25 degrees of freedom     Multiple R-squared:  0.7786,    Adjusted R-squared:  0.7697      F-statistic: 87.92 on 1 and 25 DF,  p-value: 1.157e-09      I am trying to growth of the TotalLogins (coeffient). I dont see this value in the summary. How can I get growth of TotalLogins? When I do this:               lin<-lm(TotalLogins ~ Month, data=xx)     coeff(lin)      (Intercept) Month -1673981283.6 129376.3 It is a huge number, does not seem right. Any ideas?


title: how to calculate baseline hazard rate for proportional hazard model with time varying covariates

text: I fitted a proportional hazard model with time varying covariates. However SAS does not produce the baseline hazard rates for the model with time varying covariates. I am wondering if there is an easy way to calculate baseline hazard rates so that I could use the baseline hazard estimates and coefficients to do prediction for a new covariate pattern. Your help on this is appreciated.


title: calculate Levenshtein Distance for web click stream data

text: I want to calculate Levenshtein Distance between to web click paths. I have a web page list around 500.               user id     path     1  abc,cde,eg,ba     2     abc,cde,ba     3   abc,yz,zx,eg     4  abc,cde,eg,ba     5        abc,cde     6            abc     7      cde,eg,ba      I converted this in to integer list like following.               id  path     1  1234     2   124     3  1563     4  1234     5    12     6     1     7   234      But the problem with my situation is I have 500 web pages. How to differentiate page 3 and 7 and 37? Is there more elegant way to convert this?


title: Accounting for groups in GzLM Poisson regression

text: I am currently trying to perform a regression on data with one predictor and one covariate. The hypothesis I am trying to test is whether the incidence of pores (a hole through a cell) increases when I stretch the cells. I have three cell lines and three stretch levels (0, 10 and 20%). At each stretch level and for each cell line I have 4 (and sometimes 3). Because I am sampling the occurrence of rare events (typical counts between 0 and 20) over an area, we are looking at a Poisson process. This is the reason why I assume the data is Poisson-distributed and not normally distributed. x1=1x34 vector containing the strain-level for each specimen x2=1x34 vector containing the cell line fore each specimen y=1x34 vector containing the pore counts for each specimen There are 3x3x4-2=34 specimens: cell line 1 is missing one specimen at 10%, cell line 2 is missing one specimen at 20%. I am doing my statistical analysis in MATLAB. Before I came to the conclusion that my data is Poisson distributed, I used ANCOVA with the aoctool function, fitting separate lines through each cell line. The aoctool then calculates an F-statistic and a corresponding p-value.               [h,atab,ctab,astats]=aoctool2(x1,yT,x2,0.05,'Strain','Pore Count','Cell Line','off','separate lines');      The aoctool fits the data to the following model: y = (α + αi) + (β + βi)x + ε where y is the dependent variable, α the general intercept, αi the cell line specific intercept, β the general slope, βi the cell line specific intercept and ε the error term. Fitting the data to this model regresses the pore count with the stretch-level whilst accounting for possible differences between cell lines (groups). Now I want to do the same with the fitglm (or glmfit) function in Matlab but using a Poisson rather than a normal distribution. However, I have not been able to get the fitglm function to fit separate lines for each cell line. Rather, it takes the cell line as a covariate. My question is twofold:   1. **Is there a difference between using cell line as a covariate or as grouping variable?**   2. **How do I adjust the fitglm-input to take into account cell line as a grouping variable?** This is what I have tried so far:               tbl=table(yT,x1,x2); % enter my data into a 'table' format for the fitglm function     tbl.x2=nominal(tbl.x2); % x2, cell line, is a nominal variable with values 1, 2 or 3     model=fitglm(tbl,'yT ~ (x1+x2)','Distribution','Poisson')      this results in Estimated Coefficients: Estimate SE tStat pValue   **____** **_____** **___** **______**               (Intercept)      2.5086     0.097237     25.798    9.234e-147     x1             0.024509    0.0024828     9.8714    5.5397e-23     x2_2           -0.47372      0.10328    -4.5869     4.499e-06     x2_3           -0.65656      0.10444    -6.2862    3.2532e-10      Thank you so much for your help!


title: statistical method

text: I want to compare time spent in shopping between two group people. One group has 1426 cases, whereas the other group only has 96 cases. Which statistical method should I use, t-test or Mann-Whitney test?


title: Implementing R Machine Learning Model in Real World

text: I've been able to create some successful R ML models using some of the popular machine learning algorithms. However, I'm not sure how to implement the model where the end users (technically challenged) can utilize the model with a few button clicks. They prefer a web or a traditional application interface. Can someone point me in the right direction? Thank you very much in advance.


title: Comparison of a numerical and a stochastic maximization of a Cauchy likelihood

text: I'm trying to provide a comparison of a numerical and a stochastic maximization (using uniform sampler) of a Cauchy likelihood in terms of the sample size via drawing a plot. My problem is associated with the stochastic part. For this part of R code, I don't know which function would be appropriate to find the maximum of theta to show stochastic process suitably (Of course optimize function must not be used). My R-code is as follows:               #numerical method     set.seed(1)     xm=rcauchy(500)     f=function(theta){-sum(log(1+(x-theta)^2))}     for (i in 1:500){     x=xm[1:i]     theta1=optimise(f,interval=c(-5,5),maximum=TRUE)$maximum     print(theta1)}          #uniform sampler     set.seed(1)     xm=rcauchy(500)     theta=c(runif(500,-5,5))     f=function(theta){-sum(log(1+(x-theta)^2))}     for (i in 1:500){     x=xm[1:i]     theta2=  #Related to my question     print(theta2)}          plot(theta1,theta2,xlab="theta from numerical method",     ylab="theta from uniform sampler",lwd=3)      Please feel free to edit the uniform sampler part of R-code if there are any mistakes. And I am specifically looking for a plot like this one by using the above code to make the comparison possible. ![enter image description here](http://i.stack.imgur.com/bk1ki.png)


title: Interpretation of the area under the PR curve

text: I'm currently comparing three methods and I have the Accuracy, auROC and auPR as metrics. And I have the following results : Method A - acc: 0.75, auROC: 0.75, auPR: 0.45 Method B - acc: 0.65, auROC: 0.55, auPR: 0.40 Method C - acc: 0.55, auROC: 0.70, auPR: 0.65 I have a good understanding of accuracy and auROC (to remember well i often try to come up with a sentence like "auROC = characterize the ability to predict the positive class well", while not exactly correct it helps me remember). I have never had auPR data before and while I understand how it is built I can't get the "feeling" behind it. In fact I fail to understand why the method C has an incredibly high score for auPR while being bad/average for the accuracy and auPR. If someone could help me understand it a little better with a simple explanation that would be really great. Thank you.


title: Music classification into genres

text: I need to classify music (songs) into genres (rock, french house, trash metal, etc). My idea was to extract features from the songs (bmp, zero crossing, etc) and then apply known classification algorithms. I want to start with something simple and then improve from there if needed. What features are important for songs genre classification? What are the best scalable tools to extract those?


title: Checking for multicollinearity in a negative binomial regression model

text: I'm a beginner and have this question: In a negative binomial regression analysis, is it possible to check if there is multicollinearity? I'm trying to introduce some moderating effects (e.g. interaction terms) in my model with count data as a dependent variable, and I'd like to make sure there's no multicollinearity (I've standardized the variables of interest previously). Any help will be greatly appreciated!


title: SVM with overlapping classes

text: Following my earlier post, I have successfully run linear svm on some datasets, however for other datasets, it is failing. These are the datasets where range of one class lies within that of the next class. What I did not mention in my post was that my classes have to be ordered i.e c1 should precede c2 which should precede c3. So for example linear SVM is not very accurate in classifying the following (where c1, c2 etc are class labels)               c1={3,4,5,5,7}      c2={1,2,3,4,5,6}      c3={1,4,5,7,8}     c4={3,5,6,9,10}       linear SVM gives me the following classes               c1={5,5,5,5,5,6,6}      c2={1,1,2,3,3,3,4,4,4}      c3={}     c4={7,7,8,9,10}        What I am looking to get is                c1={1,1,2}       c2={3,3,3,4,4,4,5,5,5,5,5}       c3={6,6}      c4={7,7,8,9,10}       Is there a better method to classify such datasets?


title: hypothesis test likelihood ratio

text: Let $X_{1},..,X_{36}$ be a random sample from a bernouilli distribution with parameter $p$. Suppose the sample probability is $1/3$. Using a norml approximation and likelihood ratios, we wish to test $H_{0}:p=1/2$ vs $H_{1}:p {\neq}1/2$ at the $5$% significance level. I'm not sure what the test statistic should be. Can anyone help out? Thanks


title: Bayesian hypothesis test

text: Given a i.i.d sample $X_{1},..,X_{n}$ of bernoulli random variables test 2 hypotheses $H_{0}:p=2/3$ and $H_{1}:p=1/3$. Bayesian prior is $\pi(2/3)=1/3$ and $\pi(1/3)=2/3$. Find the bayesian criterion for acceptng $H_{0}$, find the bayesian mean square error for the test and for $n=8$ compute this mean square error using normal approximation I have found the bayesian criterion for acceptance as $\sum_{i=1}^{n}x_{i}{\geq}\frac{n+1-log_{2}(\alpha^{-1}-1)}{2}$. where $\alpha$ is a value is chosen prior to the test. How do you do the other two parts? Thanks


title: Nested cross-validation for classification in MATLAB

text: I am trying to tackle a classification problem with Support Vector Machine in Matlab using SVM. Using sample codes in Bioinformatics Toolbox documentation (SVM Classification with Cross Validation) , I am able to train a SVM and find its optimal parameters. But from Inner loop overfitting in nested cross-validation and How does one appropriately apply cross-validation in the context of selecting learning parameters for support vector machines?, I will require another outer (leave-one-out) cross-validation to ensure that the final model isn't biased. My question is, so how do I implement this nested cross-validation? I understand that the inner cross-validation is to choose the optimal parameters, but what is the purpose of the outer cross-validation and how should I do it in code? Other similar information: Nested CV with shrunken centroids and SVM On Over-fitting in Model Selection and Subsequent Selection Bias in Performance Evaluation How Wrong Can We Get? A Review of Machine Learning Approaches and Error Bars


title: Regression and link function

text: Suppose we have $E (\log (Y)) = a+bx $ vs $\log (E (Y)) = a+bx $. Can $\exp (b) $ in both cases be interpreted as a geometric mean?


title: Confused about why we would use expected value instead of MLE when estimating some parameter

text: I have a conceptual confusion about the use of the expected value of a distribution. Often, we want to estimate the most likely value of something. For example, I have X= ten observations. I know X was drawn from the poisson distribution, and I know that this poisson distribution itself has a lambda that was drawn from the lognormal distribution with mu = 4, std = .05. So given my observations, and my knowledge about the underlying distributions, I want to estimate lambda. So Bayes theorem comes in handy, with posterior = lambda given X prior = lognormal distribution likelihood = likelihood function for a poisson distribution This can be solved (e.g., with rejection sampling) to discover that our Bayes estimator for lambda is ~4.28. What we did, in intuitive/oversimplified terms, is take the "center of mass" of our posterior distribution, and found it was ~4.28. My question is: if we're trying to find the most likely value of lambda, why are we taking the "center of mass" of the posterior distribution, rather than the highest point on our posterior distribution? In this example these two are extremely similar, but you could imagine changing some assumptions in this problem, and getting a posterior distribution that is (for instance) bimodal. If this were the case, the highest point of this bimodal distribution seems like it corresponds nicely to the "most likely value of lambda," whereas the center of mass of this posterior distribution could be an extremely unlikely value of lambda (e.g., a value in the valley between the two peaks of the distribution). Any help is appreciated! TL;DR: When estimating the 'most likely value' for a parameter, why do we ever want to do so by finding the 'center of mass' (expected value) from that parameter's distribution? Shouldn't we always want the 'highest point' (MLE) on that distribution?


title: What is residual standard error?

text: When running a multiple regression model in R, one of the outputs is a residual standard error of 0.0589 on 95,161 degrees of freedom. I know that the 95,161 degrees of freedom is given by the difference between the number of observations in my sample and the number of variables in my model. What is the residual standard error?


title: Statistical test for multiple comparison with dependent observations (MATLAB)

text: I want to use matlab to perform a statistical test to find if there's a significantly higher mean return in any of n different sample populations. However, I can only find post hoc (although they can be used separately the input in matlab requires that ANOVA is run first) tests for the ANOVA, which I can't apply since the assumptions do not hold (normaly distributed, equal variance, independence between samples (do they mean between observations within the samples or between the sample populations?)). What I have is different strategies and I want to find out if a significantly higher mean return can be achieved by changing a parameter (sort of). Since the strategies are evaluated using a rolling window the returns will be dependend. Also the variance most likely differs between the samples. Details: I'm evaluating a strategy using different stocks. I want to see if I can find significantly higher mean returns by using these different groups of stocks (7 groups). The strategy is evaluated over a ~16 year period using 4 month trading periods and a 2 month rolling window. Since the trading periods are overlapping the returns will most likely be dependent (the returns depend on where the trading period starts so the returns in the overlap will not be exactly equal). However, the different samples are independent of eachother as the stocks are mutually exclusive. In the second case the stocks are not mutually exclusive. Hence the (in this case) 3 groups will be slightly dependent. Otherwise the same conditions holds here. What is the most appropriate test to test for significant different means under these conditions? How can I implement it in matlab as a standalone test?


title: What do you consider a new model versus an updated model (time series)?

text: I am having some issues explaining to [non-statistician] people that it is natural to revise the parameters of a time series (ARIMA) model if you update the model with new data (add new actual values and move the forecasting point forward). I say: The model is predicting new values of the series based on recent and older values in the series. If I have new data points and wish to forecast from time t2 rather than point t1 (t2 > t1) then I have to update the model with data between t1 and t2, otherwise I am really forecasting from t1. This is natural for time series. They say: Your model coefficients have changed value, therefore you have developed a new model. You should keep the coefficients the same even with the new data. The reason for this debate is that I was instructed not to build a new model. In their opinion, I have built a new model. In my opinion, I have not. Semantics aside, I have three main questions: (1) Am I wrong? I freely admit that the model may no longer be the "best" model: whether it is the best or most suitable is not at issue. Fundamentally, the question is about when a model should constitute a "new" model. (2) Is there a practical, not-fundamentally-flawed way to do what they ask, namely to use the time series model with more recent data without changing the coefficients? The only ways I can think of are: (a) just forecast from the original point but use actual data where actual data exists, and (b) chop out the time between t1 and t2, so that the forecast period begins at t1. (I should point out that this is technically an ARIMAX model because I use a forecasted macroeconomic variable as an input.) (3) Assuming I am not wrong, does anyone have a source that might help explain time series to people who know zero statistics.


title: how to calculate E[vech(x x')vech(x x')']?

text: Supposing a vector x follows normal distribution. I want to calculate the expectation of the "fourth moment" in a vector form, meaning $\text{E}[\text{vech}(x x')\text{vech}(x x')']$, given that we know $\text{E}[x x'] = S$ and the expectation of $x$ is a zero vector. This is a piece of a Fisher information matrix. Anyone can give me some suggestion on how to get that? Thanks for your help.


title: Best estimate for a decile or quintile mean from a known distributional family?

text: Suppose you have a population drawn from a known distributional family f with a vector of unknown parameters θ. You may also assume that the distribution is strictly non-negative and skewed, with a finite mean and variance. Suppose further that you have a large sample, sorted least to greatest. The sample is then divided into n groups by n-tile (e.g. 10 groups by decile). For convenience, assume that the number of observations is divisible by n. You use the entire sample to produce a MLE of the unknown θ. From this point, one could estimate the quantiles numerically by finding the values for which the PDF f(θhat) is equal to an integer multiple of 1/n, and then the quantile group means by integrating xf(θhat; x) between adjacent quantiles. Or, you could find the quantiles by dividing the sorted population into n groups of equal size and taking the sample quantile group means as estimates of the true quantile group means. If my interest is in the quantile group means, rather than the quantiles themselves, is there any good reason to prefer one of these estimates over the other, in terms of bias, variance, MSE, or some other standard basis? Or does it depend on the underlying distribution, and if so, how? If anyone wants to talk about any other candidate estimator, I'd be happy to hear that as well. Any help or pointers to literature will be greatly appreciated.


title: Survival analysis, one cohort, two classifications

text: I'm wondering about this question: short version: how to adequately compare the effect of a reclassification of the same subjects on survival long version: I have one cancer cohort that was sorted into TNM classes long time ago. The definition of the classes has now been updated and I like to know if, say, class 1 of the old classification has a better/worse survival than class 1 of the updated classification. My problem is that the cohort stays the same, only the sorting into classes changes. Mini example: patients 1,2 are in class I_old and patients 3,4 in class II_old, patient 1 stays in class I_new but patient 2 goes to class II_new, patient 3 stays in II_new, patient 4 changes to I_new. If I want to compare the effect of the reclassification in class I, I would have to compare the survival of patients 1 and 2 versus the the survival of patients 1 and 4. So patient 1 is in both of the groups. I doubt that the standard procedures like logrank test are adequate here. But what can be done? I have a very vague idea about proportional hazards regression with a time varying covariate (namely indicator for "before" and "after" update) or a frailty model but cannot tell at all if I'm going in the right direction. What do you think? psj


title: Statistics and data mining software tools for dealing with large datasets

text: Currently I have to analyze approximately 20M records and build prediction models. So far I have tried out Statistica, SPSS, RapidMiner and R. Among these Statistica seems to be most suitable to deal with data mining and RapidMiner user interface is also very handy but it seems that Statistica, RapidMiner and SPSS are suitable just for smaller datasets. Could anyone recommend a good tool for large datasets? Thanks you!


title: Model suggestion

text: Could anyone give me hints as to a model framework that can be used in the following setting: The outcome A is dichotomous. I want to investigate the effect of a continuous variable B and a continous variable C on A in a longitudinal setting. So far so good. The problem is that C depends on B, in that values of B above a threshold (Unfortunately we do not know the threshold and the threshold itself might not be constant over time) will change C. The change in C will furthermore reduce B in the next measurement. How can I decipher the isolated effect og B and C on A? Would a mixed-logistic-regression model be OK? ie;t=time A~B*t+C*t+B*C*t Any other ideas? //M


title: Proportion data: - beta regression or logit transformed OLS regression?

text: I have proportion data (percentage viewership of TV programs) that i'd like to model as a function of various demographics (age, sex etc.) and time (year). After surveying options for appropriate multiple regression models, I'm debating between the following two strategies: 1) fit a beta regression model after dividing the percentage data by 100 and adjusting the range slightly so that values of zero and one do not occur. 2) fit an OLS model after logit transforming the percentage data (again, divided by 100 and adjusted slightly) so that the dependent variable is mapped to the Real line. One key consideration is that i'd like to make the results as intuitive as possible to a non-statistical audience. So, interpretations such as "for a one unit change in X we get a percent change in Y", or something like that, would be most welcome. Can anyone outline pros and cons of these two approaches in this regard? It seems to me that using beta regression with a logit link, and then calculating odds ratios may lead to nice "percent change" explanations. The coefs from the OLS model would also be on the ln(odds) scale, so I assume I could also do the same for that model. My data are too large to share, but I ran both models in R and there are only minor differences in the coefs.


title: How to get the confidence interval around the variable importance generated by randomForest?

text: Using randomForest in R is it possible to get a variance or confidence interval around the importance of a variable (% difference in mean square error)? From the randomForest package               set.seed(4543)     data(mtcars)     mtcars.rf <- randomForest(mpg ~ ., data=mtcars, ntree=1000, keep.forest=FALSE,                               importance=TRUE)          importance(mtcars.rf, type=1)      > >     %IncMSE >     cyl  17.061579 >     disp 19.020929 >     hp   18.195448 >     drat  6.677707 >     wt   17.858616 >     qsec  5.849071 >     vs    5.467722 >     am    3.933137 >     gear  5.129811 >     carb  8.332738 >                 varImpPlot(mtcars.rf)      ![Variable Importance Plots](http://i.stack.imgur.com/WyvLk.jpg)


title: What's wrong with Johnson SU fit using moment=find (SuppDists)

text: I have a vector of returns `r` that have skew and kurtosis. I want to fit a Johnson SU distribution to this data. Using the code below:               plot(qJohnson(ppoints(r),JohnsonFit(r, moment='find')),     sort(r), xlim=range(-.2,.2), ylim=range(-.2,.2)); abline(0,1)      This results in the QQPlot below. Clearly there is something wrong with the estimation. When I use the `moment='quant'` I get a better fit (most of the points are over the x=y line. Am I using the JohnsonFit function right? ![JSU Fit](http://i.stack.imgur.com/grSTn.png)


title: Marketing research survey results

text: I added a question about marketing research survey result. There are two questions and two results related to this survey. Can you help me answer these questions about this data sheet. ![http://i.stack.imgur.com/9NQYR.jpg](http://i.stack.imgur.com/jn4gN.jpg)


title: Confidence interval definitions

text: Wikipedia provides the following definition for a confidence interval for a parameter $\theta$: > **A confidence interval** for the parameter θ, with confidence level or > confidence coefficient γ, is an interval with random endpoints (u(X), v(X)), > determined by the pair of random variables u(X) and v(X), with the property: > > ${\Pr}_{\theta,\varphi}(u(X)<\theta<v(X))=\gamma\text{ for all > > }(\theta,\varphi). $ > > Here $Pr(θ,φ)$ indicates the probability distribution of X characterised by > $(θ, φ)$. > > In a specific situation, when x is the outcome of the sample X, **the > interval $(u(x), v(x))$ is also referred to as a confidence interval** for > $θ$. Note that it is no longer possible to say that the (observed) interval > $(u(x), v(x))$ has probability γ to contain the parameter $θ$. This observed > interval is just one realization of all possible intervals for which the > probability statement holds. I presume that $X$ here is the random variable representing the sample and that $u(x)$ and $v(X)$ refer to not only random variables, but to **methods** for constructing the upper and lower bounds of the confidence interval from one sample. So it looks like there are **two** confidence interval definitions. In the first one, one can say that the interval has probability $\gamma$ of containing $\theta$, while in the second one, one cannot. Why is this the case? Perhaps more specifically, I fail to see how one reaches the following conclusion **given the second definition of a confidence interval** provided: > ... Note that it is **no** longer possible to say that the (observed) > interval $(u(x), v(x))$ has probability $\gamma$ to contain the parameter > $θ$.


title: What's the difference between "deep learning" and multilevel/hierarchical modeling?

text: Is "deep learning" just another term for multilevel/hierarchical modeling? I'm much more familiar with the latter than the former, but from what I can tell, the primary difference is not in their definition, but how they are used and evaluated within their application domain. It looks like the number of nodes in a typical "deep learning" application is larger and uses a generic hierarchical form, whereas applications of multilevel modeling typically uses a hierarchical relationships that mimic the generative process being modeled. Using a generic hierarchy in an applied statistics (hierarchical modeling) domain would be regarded as an "incorrect" model of the phenomena, whereas modeling a domain-specific hierarchy might be regarded as subverting the objective of making a generic deep learning learning machine. Are these two things really the same machinery under two different names, used in two different ways?


title: Computing the Correlation with Skewed Repeated Measures Data

text: Suppose I have two variables (e.g., Blood Pressure and Cholesterol) and I want to investigate their association across time. Also suppose that each individual has been measured three times (longitudinally) on each variable. So, the data might look something like this:                 ID Time  BP  CH     1  1    1 124 214     2  1    2 140 218     3  1    3 135 209     4  2    1 126 206     5  2    2 132 197     6  2    3 128 204      My collaborators (non-statisticians) just computed a spearman correlation between the two (ignoring the dependence). They used a Spearman because the data are highly skewed. This is obviously a problems since the data are not independent. I saw this article that shows how to compute this using mixed modeling. However, my distributions are pretty messed up, so I was thinking about combining the approach mentioned earlier with this approach (which summarizes a nonparametric extension of the mixed model). So, my questions: (1) Am I reinventing the wheel? Has someone already tackled this problem (and if so, is there an R package or something)? (2) If this hasn't been done, does this seem sensible (first estimate a non-parametric mixed model, then use the estimates from that to compute the correlation between the variables)?


title: Linear-logarithmic Regression in MATLAB with two input arguments - which model to assume?

text: Look at the following plot and ignore the solid lines please (just look at the dotted/dashed ones). ![enter image description here](http://i.stack.imgur.com/0QaAa.jpg) ONE curve is described by the following parameters and ONE fixed value for the second input parameter V:   * X-axis: Input-Parameter `T`   * Y-axis: Output-Parameter `g` For each curve, `g` is between [0, 255] (thus always positive), concave, bijective. I know from the process that lies behind the measures, that by increasing `V`, the corresponding curve flattens. The different curves result when varying `V`. The orange curve at the top is for like `V=100`, the bottom curve (red/magenta) results for `V=180`. I have measured data with a lot more data points in the following form:               T[1]  V[1]  g[1]     T[2]  V[1]  g[2]     T[3]  V[1]  g[3]     ...   V[1]  g[4]     T[N]  V[1]  g[5]     .......     T[1]  V[N]  g[1]     T[2]  V[N]  g[2]     T[3]  V[N]  g[3]     ...   V[N]  g[4]     T[N]  V[N]  g[5]      Now I want a regression like this:               g = g(V, T)      which would yield the curve for a fixed V-value:               g = g(T), V=Vfix      Which regression-funktion in MATLAB do you think would work out the best way? And how to assume a "model" here? I only know (from the process itself AND obviously from the plots), that its some sort of linear curve at the beginning, pass over into a logarithmic curve, but I dont know how the value of `V` inferfers with it!? Thanks a lot in advance: for any advice..


title: Compute how much the one percent wealthier concentrates with negative data

text: I want to compute, using the Survey of Consumer Finances (SCF) database, how much the one percent wealthier concentrates. The variable I'm using to measure wealth is the "networth". The problem is that this variable has negative values and, as far as I know, we cannot do this calculation with negative values. The result could be negative, which makes no sense. Anyone knows how could I get around this? Is there any transformation that I can make in the data that passes over the problem without changing the essence of the distribution? Thanks!


title: The probability density function of a rescaled / transformed chi-squared random variable

text: If $X\sim\chi^2_{6}$, what is the probability density function of $T = \frac{(X-6)}{\sqrt12}$? The problem I'm confronted with is that the chi-squared random variable, $X$, can assume only positive values while the random variable $T$ can assume both positive and negative values. I tried to use both the PDF technique and the CDF technique for transformations of random variables but in the end the PDF of $T$ could assume only positive values, that is my problem.


title: Predicting the impact point of a moving object

text: Suppose we have a moving object (a horizontal projectile motion as one of the most basic examples). Is there any way to predict where it will hit finally? Please note that I'm looking for a machine learning method not a closed form solution. Although we can track the motion, using Kalman filter, that is only applicable when we want to predict the new future _(As far as I'm considered)_. But I need to predict the ultimate goal of a moving object. To better express the problem let see the following example: Suppose a goalkeeper robot that of course uses filtering methods to smooth the ball motion. It needs to predict if the ball is going to enter the goal or not, before it decide to catch the ball or neglect it to go out. Thanks for any help


title: Determining number of tests to correct for using Bonferroni

text: I am doing final revisions to my dissertation, which includes a large number of regressions. I have six outcome variables (DVs) and I am testing the effect of four primary (correlated) IVs, but I am testing them separately to determine the unique contribution to the model of each IV. Therefore, I have 24 regressions. One of my committee members would like me to use Bonferroni correction to control for Type I error; however, I am struggling with determining the number to use to determine the significance threshold. I guess the main thing I'm struggling with is that while I include these 24 regressions in my dissertation, any publication that is derived from my dissertation will not include all of these analyses, so setting the Bonferroni correction to .05/24=.002 seems (1) extremely conservative and (2) would be less conservative in a subsequent journal publication, as I would be correcting for fewer tests. Are there other ways to approach determining the number of "tests" that alleviate Type I error concerns and make theoretical sense without being so conservative? I've been struggling to find some useful resources/articles for the last few days that account for a situation like mine, but most discussions of Bonferroni (or Holm-Bonferroni or other derivatives) focus on treatment of a single regression or analysis rather than a range of analyses across a dataset. I also realize that this is one of the biggest criticisms of Bonferroni--that the idea of constantly adjusting the alpha with each test would (in theory) require us to put off publishing anything until we had exhausted a given dataset. Any links to useful articles or advice on how to approach my current problem would be greatly appreciated!


title: Regression: Finding the vector explaining most variance in 3D space

text: I have a specific regression problem of which I am not sure how to solve it. We have a dependent variable which varies depending on a certain area that is electrically stimulated in the brain. The stimulation location is represented as 3D coordinates. A multiple regression with the three coordinates (i.e. axes) as predictors yield that two axes are significantly correlated with the DV. Hence, we should be able to find a vector that explains most of the variance and is not one of the „c ardinal“ axes but somehow diagonal. Can anybody recommend any statistical methods (maybe somehow similar to PCA) to obtain this vector? I am searching for an automatic way of finding the direction that is decisive in modulating (i.e. explains most of the variance) the DV. Thanks for any advice! Petra


title: Creating data that follows a specific data distribution

text: I have a variable z which has around 3000 values between 0.0 and 0.5. I have plotted here z in the y-axis and 3000 evenly spaced numbers over the interval 0.0 and 0.5 in the x-axis so that one can see my data distribution. ![enter image description here](http://i.stack.imgur.com/d9DMf.png) I want to create another set of data of 100,000 points, **that follows the distribution of my original data, i.e z**. How can I do this? Preferably in Python?


title: Multinomial Naive Bayes performance with and without Bigrams/Stopwords

text: I have a dataset of 300 text documents that I have manually classified into three classes. The first two classes are relatively equal in size; the third class is about 40% larger than each of the individual first two classes. The classes are very similar, i.e. they share similar features. After pre- processing and removing words that only appear in 2% of the texts, Multinomial Naïve Bayes predicts about 48-58% of the test set (20% of the texts) correctly. I notice that the larger third class is predicted more accurately than the smaller first and second classes. I would like to refine pre-processing to achieve 65-75% accuracy. The first thing I did was to tokenize the features into bigrams. After doing so, the model only predicts between 40-50%, instead of 48-58% with single words. Intuitively, I expected the model to perform better, not worse, with bigrams. Shouldn't bigram tokenization reduce the amount of shared tokens across the classes and increase the predictive power of the non-shared tokens as a result (at least marginally)? I am a bit concerned because my next refinement is to create a stopword list of words that frequently appear across all classes but have no predictive content. My reasoning was that by reducing the amount of high-frequency shared words, I would increase the predictive power of the remaining words. Isn't this the point of stopword lists?


title: Meaning of output terms in gbm package?

text: I am using gbm package for classification. As expected, the results is good. But I am trying to understand the output of the classifier. There are five terms in output.               `Iter   TrainDeviance   ValidDeviance   StepSize   Improve`      Could anyone explain the meaning of each term, especially the meaning of **Improve**.


title: Association Rules "with a kind of class"

text: I want to use/adapt a recommendation algorithm for posters in an e-commerce. The thing is that I want to use previous categories searched before posting in a particular category (has to be at a very granular level), so I need to create rules that always comes up with the posting/replying class. So I have 2 different kinds of item: viewed and posted/replied. I was thinking to create rows with **previous browsing data** (maybe 48hrs before posting) and then a variable dummy for each category of posting/replying. But don't know if this will be the correct way. Also I would like to consider previous posting/replying and in that case I will need to consider a more time (2 or 3 **previous months posting/replying data** ). Do you have any ideas of how to prepare the data to run a recommendation algorithm for this case? Thanks!


title: Does it make sense to fit a Pareto Curve to sales data?

text: The Pareto principle is used surprisingly widely in business. I'm wondering how correct it's wide use is. It seems as though it is often used without empirical verification as if the situation fits a Pareto distribution a priori. Here's an example: Let's say there is a retail store with a wide range of items in a single aisle. The items are assumed to be a Pareto distribution. The sales data of the items are then fit to a Pareto curve and shown to fit it (using R^2, a basic measurement). This curve would in turn help you find out, I suppose, what the best items would be to take out and so on. My question is two-fold: 1) Is it actually so easy to empirically verify whether the aisle's historical sales is a Pareto distribution? 2) Does fitting the data to a Pareto curve mean anything? I've seen that, broadly, fitting data to a Pareto curve does not necessarily imply that the data is actually Pareto distributed.


title: LIBSVM RBF kernel time complexity

text: Is it known what is the LIBSVM time complexity when using RBF kernel functions?


title: Measure precision

text: This is used in information retrieval. I have an algorithm that uses a sample set to predict a yes / no. I think the correct term is binomial. Yes this document is about sports or no this document is not about sports. I need to represent the effectiveness of the algorithm with that population and that question.   Can I just compare the output from algorithm to the sample set. Clearly I need remove that item from the sample set to fairly evaluate the prediction. If I have a sample set of 1000 and as reviewed by an expert exactly 700 were about sports. If I go through the sample set one at a time and remove that item from the sample set and predict and compare the prediction to the actual value is that a fair measure of precision. If the prediction was correct on 990 of the items can I say the precision is 99%? I think it is silly to have have to get another X samples to test the prediction. If I remove that item from the sample set then can I not pretend it was the test of the sample set / algorithm? The sample set is random. I also use the sample set to estimate the number in the total population. I get that if I tweak the algorithm for conformance to the sample set then I cannot use that as a measure of precision. At that point I would need to test against a blind sample.


title: What are some linearly inseparable data sets for testing support vector machines and artificial neural networks?

text: I am looking for some linearly inseparable classification problems and nonlinear regression problems. What are some public data sets?


title: Poisson GLM vs Quasi Poisson GLM

text: I have a Poisson GLM in R that is over dispersed, so I fit a quasipoisson GLM, however the residual deviance nor the degree of freedom change. Can that happen? What does it mean in that case? Thank you.


title: Level of analysis and assumption of independence in Pearson Chi-square test

text: I have survey data about injuries occurring to athletes over the span of their careers. Many of the athletes suffered multiple injuries (the most any athlete reported was five and this did not occur very often). While I can see that injuries within an athlete may be dependent on a previous injury, in the data, many within-athlete injuries involved different body parts with different mechanisms and types, occurring years apart. My interest is in the characteristics of the injury itself not the athlete. I am looking to see if any one particular mechanism of injury was prevalent in the sport. And with that information, look at possible modifications to the equipment or rules to improve the safety. I used a Pearson Chi square test to see if Anatomical site of injury was distributed evenly across mechanism of injury (e.g. contact with opponent, quick turn, etc.); and the type of injury (e.g., strain, sprain, etc.) across mechanism, too. My original intent for this data was to present descriptives and some graphs contrasting the different mechanisms and the parts injured or the types of injury, however, I was asked by a reviewer to perform some statistical analyses to show that, for example, shoulder strains occurred more frequently than shoulder lacerations. Now I have a concern about the independence of the observations. Since I'm looking at injury level variables and each injury is only included once in each analysis, I don't think I'm violating that assumption despite several athletes being in the analysis multiple times. Is this correct? Or do I have to worry about independence. There are a number of pieces of equipment that can be modified and this will be the next project once "high risk" equipment is identified. With this project I was hoping to get evidence that would justify targeting those pieces of equipment with some data that reflect higher injury rates when using specific equipment and thus avoid altering equipment that did not have a high incidence of injury associated with it. Thoughts and advice appreciated.


title: Scores derived from variance specific for a factor

text: I received a comment about my article from one of the reviewer but I do not understand what I should do. Here is a short explanation of the study: This study investigated the influence of lectures about the Renaissance and abstract art on ratings of paintings from these two periods in art history. The study included two sessions. In the first, 72 naive participants rated on 21 seven-point Likert type scales the figural and abstract paintings. In the second session participants were divided into three groups and they rated a new, parallel set of paintings (on the same scales) after attending a lecture on Renaissance or abstract art or without lecture attendance. Three dimensions were extracted: Aesthetic experience, Relaxation tone and Arousal. On these three dimensions I explored effects of different experimental variations (different lectures). Three-way mixed ANOVA was used with 2 within-subjects factors: Time (2 levels, Pre- and Post-treatment) and Style (2 levels, figural and abstract) and one between-subjects factor Lecture (3 levels, lecture on Renaissance paintings, abstract paintings and control group without lecture). I received following comments: Since comparing pre- and post-lecture scores per factor showed generalized effects across factors, I am wondering whether these effects originate from the (substantial) shared variance of factors, i.e., from the G factor of subjective experience. It would be interesting to test effects of lectures on scores derived from variance specific for a factor. Should the lecture effects disappear one might argue that lectures influence the G rather than each factor separately. I am not sure how to get perform what is suggested. Thanks in advance


title: How to correctly concatenate time series data

text: I have this shifted time series data. For one set consists of features from week1-5 and labels at week6. Another set features from week2-6 and labels and week7 and so on. I have like four sets of them. If I consider the individual set, then the positive labels are very scarce. So I was thinking of concatenating these four sets so that I have abundant labels. However, in each dataset, the features are almost the same except for a week. So lets say this particular example, lets say 5th in each of the four datasets which belong to the same user has positive label in each of the four datasets We can say that this example, is repeated four times or weighted four times if I concatenate all the four sets. This may lead to overfitting as well. What is the best way to handle this issue?


title: Error on extrapolated values from a fitted function

text: I have some data points and I fitted a function (2nd order polynomial here) to the data. The algorithm (`scipy.optimize.curve_fit`) gave me the optimal parameters and the covariance matrix of those parameters. In the documentation, it says that the square roots of the diagonal entries give the errors for the parameters. Now I that I have the function $f$ fitted to points $0 < x$, I would like to find the value of $f(0)$ and its error. $f(0)$ is easy, but how do I acquire an estimate for the error?


title: Can I hypothesise an insignificant mediation? (Very basic stats)

text: This is for a simple psychology report. I was wondering if I could hypothesise an insignificant mediation. Eg: "...based on this research it is predicted that group identity is not a significant mediator (or does not sognificantly control) narcissism's relationship with aggression." Then I would run a mediation analysis using indirect syntax for SPSS and report that narcissism has a significant relationship with aggression, which is not significantly influenced by (does not significantly overlap with) group identity. Or should i just include it in a regression and that would be clearer? But then my main question remains, can I predict in my hypothesis that a variable will have an insignificant correlation value? Or a weak relationship? I basically want to say that i predict that narcissism will interact with rumination on aggression (using modprobe syntax and correlational data) . But i would also like to add in my initial prediction (based on research) that narcissism's relationship with aggression has nothing to do with group identity. Thank you all for any input. As I said this is a basic psychology report, not worrying too much about null hypothesis and too much other hardcore statistical stuff. I just need a basically worded prediction (which i think i have) and a basic way to test it. So the basic question is: can I speciffically predict insignificant correlations, mediations or even interactions? And would confirming this prediction support the notion that the two are unrelated in their relationship with aggression? (Or even that group identity has no relationship with aggression, not even via narcissism?) Thank you for reading.


title: Studentized residual distribution

text: I read that in a regression with $k$ regressors, the t-statistic corresponding to a certain coefficient follows a $t(n-k)$ distribution. However, later on I read that studentized residuals follow a $t(n-k-1)$ distribution. How could it be that the extra degree of freedom is lost: isn't this also just based on a regular t-statistic?


title: Cross-Section Dummy Variables in a Logistic Regression

text: first of all - I am not really into regression analysis (I have performed some linear regressions) and my statistic/econometric skills are more or less elementary. I have panel data (several countries over a certain time period) and I would like to apply a formula such as this: ![enter image description here](http://i.stack.imgur.com/XUFH3.png) Source: Button, Ngoe and Hine (1993) The authors state that this is a "quasi"-logistic function. Now is there some sort of literature on logistic/quasi-logistic models that is easily understandable? Especially: how do I calculate these country specific constants c? All that is said in the paper on this part is: > Cross-sectionally, countries vary in many ways which are impossible to > quantify. In order to allow for these local features a set of dummy > variables (c), one for each country, is included in the model. These take > the values of unity if an observation relates to a country and zero > otherwise. They are specified so that th error term varies across countries. Thank you very much!


title: PCA-based dimensionality reduction when the number of data points is less than the number of variables

text: Assume we are given a $p \times n$ (variables $\times$ data points) data matrix $X$, with $p > n$ (i.e. more variables than data points). Performing PCA on such a matrix yields an $n \times n$ eigenvector matrix $U$ (instead of the desired $p \times p$, because $p > n$). Now suppose we use the eigenvalues vector $\lambda$ associated with $U$ to choose the top $k$ eigenvectors that retain, say, $99\%$ of the variance. We then set all other $n-k$ dimensions to zero, obtaining a matrix of the same dimensionality with $U$ ($n \times n$), that contains zero columns. Let's call this matrix $U_{reduced}$. Normally, if $p < n$, to reduce the dimension $p$ of $X$ to $k$ (in the principal components space), we would do the following computation: $X_{rotated} = U_{reduced}^T X$, and we would retain the top $k$ (non-zero) rows of $X_{rotated}$. Now that $p > n$ the dimensions in the above computation do not match: $(n \times n)(p \times n)$. To "resolve" this I increased the rows of $U_{reduced}$ from $n$ to $p$ by padding them with zeros. Is this the correct way of doing this? Thank you.


title: Convergence in distribution, probability, and 2nd mean

text: > Let $\mathbb P(X=1) = \mathbb P(X=-1) = 1/2$. Define > > $$X_n = \begin {cases} X & \text{with probability } 1- \frac{1}{n}\\\ e^n & > \text{with probability } \frac{1}{n} \end {cases}$$ > > Does $X_n \leadsto X$ (i.e., converge in distribution)? Does $X_n > \stackrel{P}{\rightarrow} X$? Does $\mathbb{E}[(X-X_n)^2] \to 0$? I think $X_n \leadsto X$ because just by picturing the CDF as n goes to infinity, their CDFs converge. Can I get some confirmation on this? But for convergence in probability, I'm having trouble with the limits. To do this, what I'm trying to do is find $\mathbb P(|X_n-X| \leq \epsilon)$ and then just take 1 minus that as shown on this site: http://www.statlect.com/prbcon1.htm If $X_n = X$, then $|X_n - X|$ is just $0$ so for all $n$, $|X_n -X| < \epsilon$ But if $X_n = e^n$ then I have the two cases for $X$. For $X=1$, $|X_n-X| = |e^n-1| \leq \epsilon$. and For $X=-1$, $|X_n-X| = |e^n +1| = e^n + 1 \leq \epsilon$. This is where I am stuck. I'm trying to find conditions relating $n$ and $\epsilon$ similiar to here the example on the site. But I'm not sure how to get that from $|e^n - 1| \leq \epsilon$ and $e^n + 1 \leq \epsilon$ . Am I even on the right track?


title: How to account for multiple & varying amounts of observations per factor level & still retain info? As a random effect in GLMM or take the mean?

text: I would appreciate any help! Specifically I would like to know which option is best. _Question_ Does Var2 influence Var1 in relation to the factor(Ind), and does Var3 & Var4 also have some effect. **The problem** I have multiple observations of Var1 per factor level (or per Individual) and the number of total observations (the counts), varies a fair bit per individual. But with my other variables 2-4 I have only one observation per individual. **Option 1:** Here I just repeat Var2-4 and can take every observation of Var1 into account.               head(df1)       Ind  Var1 Var2 Var3 Var4     1   A   2.8  1.3  4.3    F     2   A  14.2  1.3  4.3    F     3   A  10.4  1.3  4.3    F     4   B  22.1 -0.1  2.0    M     5   B 135.5 -0.1  2.0    M      So to account for this I put Individual as a random effect in my mixed model. mod1 = lmer(Var1 ~ Var2 + Var3 + Var4 + (1| Ind), data = df1) My predicted variable (Var1) is normally distributed once transformed. **Option 2:** To make sure I don't lose information, but have just one row per individual, I take the mean, min, max and counts of the observations. Please recommend any measure you find more appropriate.               head(df2)       Ind Var1_Mean Var1_Min Var1_Max Var1_Count Var2 Var3 Var4     1   A      9.10      2.8     14.2          3  1.3  4.3    F     2   B     15.11      9.2     20.9         17 -0.1  2.0    M     3   C     37.70      2.9    133.2         33  1.5  3.7    F     4   D     38.10     22.8     61.1          4  1.2  2.9    M     5   E     39.40      5.3     87.5         15  0.9  2.9    M      mod2 = lm(Var1_Mean ~ Var2 + Var3 + Var4, data = df2) I do the same as above for Min & Max and all are transformed so they are normally distributed. BUT for Var1_Count mod3 = glm(Var1_Count ~ Var2 + Var3 + Var4, data = df2, family = poisson) I have run all models, and there are differences, generally not so great. But I want to make sure I am addressing this in the right way. Which option do you think (if any) is best? Thanks in advance! To note: I am new to both GLMMs and R. But keen to learn.


title: Relation between statistical randomness, uniform distribution and independence

text: In Monte Carlo simulation, we often consider how well a sequence of generated points are. If I am correct, one aspect is statistical randomness: > A numeric sequence is said to be statistically random when it contains no > recognizable patterns or regularities; sequences such as the results of an > ideal dice roll, or the digits of π exhibit statistical randomness. I was wondering if a sequence is statistically random, then does it mean that   * when the points in the sequence are viewed as realizations of a random variable, the random variable has a uniform distribution?   * when different points in the sequence are viewed as realizations of different random variables, the random variables are independent or even i.i.d.? Does statistical randomness have anything to do with uniform distribution and independence? Thanks and regards!


title: Difficulty in understanding statistics, inference and estimation

text: I have a random variable, X, with cdf $F$ and pdf $f$. I want to estimate a parameter of $F$, say mean, $\mu$. So what do I do? I construct an _estimator_ , $Y_n$, with several random variable X1, X2,..., Xn from $F$ such that, $Y_n=g(X1, X2, ..., Xn)$ and make sure that the random variable $Y_n$ has the property, $E[Y_n] = \mu$. Now what? Initially I had problem 1 and now I am stuck with problem 2. How will I find $E[Y_n]$? # Example: let $Y_n = 1/n*\sum_{i} X_i$ and we do have $E[Y_n] = \mu$. How do we find $E[Y_n]$? Suppose we take _single observations_ {x1, x2, ... xn} from our iid $X_i$. How will $y_n = 1/n*\sum_{i} x_i$ approximate to $\mu$? I know $\lim_{n \to \infty} Y_n$ converges in probability to $\mu$. But that does not mean the value $\lim_{n \to \infty} y_n = \mu$. Where am I wrong?


title: The role of the bias terms in matrix factorization formulas?

text: I'm reading about matrix factorization for recommender systems. A basic matrix factorization model would be something like: $(p_i \times q_j ) + b_i + b_j$. That formula would compute the rating for one element in the rating matrix. Now I can see that this is very similar to a usual logistic regression formula where you multiply the feature vector with a weight vector and add a bias terms. That would give you the formula of a hyperplane. But with the matrix factorization formula you have another bias terms. So can you explain to me what it does?


title: Result of the summary of the proportional hazard model with frailty in R

text: I am trying to fit coxph a model with frailty with syntax below:               > kidneyExt.cox <- coxph(Surv(time,inf) ~ age+sex+frailty(patient)+strata(type), data=kidneyExt)     > kidneyExt.cox     Call:     coxph(formula = Surv(time, inf) ~ age + sex + frailty(patient) +          strata(type), data = kidneyExt)                           coef     se(coef) se2   Chisq DF p           age               0.00743 0.012    0.012  0.38 1  5.4e-01     sex              -1.48647 0.371    0.371 16.01 1  6.3e-05     frailty(patient)                          0.00 0  9.3e-01          Iterations: 6 outer, 37 Newton-Raphson          Variance of random effect= 5e-07   I-likelihood = -112      Degrees of freedom for terms= 1 1 0      Likelihood ratio test=15  on 2 df, p=0.000553  n= 76       And I am curious about last line about Likelihood ratio test. Is it a result for a likelihood ratio test? And it's p-value is 0.000553? I am asking because when I type syntax like this below I got something different from summary.               > kidneyExt.cox$loglik     [1] -119.4828 -111.9835      Is that test in a summary has it's name if it's not likelihood ratio test? Thanks for any help!


title: How to derive the gradient formula for the Maximum Likelihood in RBM?

text: I am learning RBM (restricted Boltzmann machine) for deep learning. The log- likelihood of RBM is given as : ![enter image description here](http://i.stack.imgur.com/1AOEC.png) and its gradient w.r.t. the parameter is: ![enter image description here](http://i.stack.imgur.com/5smlF.png) I don't understand how is the gradient derived from the log-likelihood.


title: I do not have hypotheses for my research and I have got only categorical data. what kind of data analysis can i conduct? also

text: I do not have hypotheses for my research and I have got only categorical data. what kind of data analysis can i conduct for correlation between variables? also is chi-square used only for hypotheses testing? can it be used for exploratory analysis? and how to test if there is moderation relationship between variables?


title: Classifiers that support weighting of the instances

text: I'm thinking of implementing my own boosting algorithm, so I'm looking for any multiclass classfication algorithms that would support weighting of the examples, i.e. you could specify what examples are more important to affect algorithm optimisation. Of course I know about decision trees, but I'm looking for something else, maybe some other rule learning techniques. I want to do this in R, so it would be good if algorithms would be implemented there.


title: Generating a dataset from mean, standard deviation, n and 95% CI

text: I have the output of a couple of Socprog models and I'd like to see if the results are statistically significant. Group A output: `mean=55.8 SD=25.6 (95% C.I. 38.5–113.4), n=66`   Group B output: `mean=28.9 SD=11.8 (95% C.I. 6.0–36.2), n=116` I thought perhaps generating the rnorm dataset but that won't allow you to input the 95% intervals.


title: Need help doing multiple linear regression with several independent variables that are of differing levels of measurement

text: I am attempting to predict levels of body dissatisfaction from a number of independent variables. The dependant variable is Body dissatisfaction = overall score on a body dissatisfaction scale. The independent variables (predictors):   1. Self esteem levels = overall score on Rosenberg Self Esteem Scale   2. Drive for muscularity levels = overall score on Drive for Muscularity scale   3. Age   4. Current involvement in sport and exercise = yes/no answer   5. Previous involvement in sport and exercise = yes/no answer   6. Masculine gender role identification = overall score on Bem Sex Role Inventory scale   7. Fitness magazine exposure = amount of fitness magazines read + time spent reading these magazines.    8. Perfectionism levels = overall score on frost multidimensional perfectionism scale.    9. Levels of appearance conversations regarding muscularity = overall score on Male Body talk scale.  So basically I need to know:   * Can I do multiple regression when I have different IV's with different levels of measurement? If not, what other type of test could I use?   * Do I need to transform the nominal variables (current and previous sport and exercise involvement) in dummy variables even though there is only two possible answers (yes = 1, no = 2)?   * Can the overall score of scale items be considered interval data? If not, how can I alter it e.g. Rasch modelling I can provide more information as needed, thanks for any help in advance!


title: How do I investigate how long it takes one variable to affect another in a time series

text: I am a total newbie when it comes to time series, so it is quite possible this question is duplicated somewhere else, only that I cannot find it because I don't know what this feature is called. ### My data: I have weekly measurements of a variable "change in rank". I also have weekly measurements of another variable loosely described as "change in goodness". I believe that "change in rank" is dependent on "change in goodness", but I also believe that whenever "goodness" improves, it takes some time before rank improves. A simplistic representation of the data would look like this:               change_in_rank = c(0,0,0,0,-10,0,0,0)     change_in_goodness = c(0,2,0,0,0,0,0,0)      where this constructed example shows that a improvement of 2 in goodness took 3 time steps to lead to an improvement in rank. ### My question: What analysis method do I use to detect whether such a relationship actually exists between "change in goodness" and "change in rank", and if such a relationship exists, how do I find out how long it takes for a change in "goodness" to flow on to a change in rank. ### The motivation behind this question: Some additional, non necessary background: it is claimed that search engine optimisation takes time to affect your site's ranking in Google. E.g.: "Why does SEO take so long?" I am trying to find a statistical way of calculating how long it takes for a change in a site to actually reflect in the site's Google ranking. ### More motivation: Slightly different type of data but similar question... "How long does it take before a change in diet starts to take effect?"


title: Similarity measures between curves?

text: I would like to compute the measure of similarity between two ordered sets of points---the ones under **User** compared with the ones under **Teacher** : ![Example](http://i50.tinypic.com/2uj0580.jpg) The points are curves in 3D space, but I was thinking that the problem is simplified if I plotted them in 2 dimensions like in the picture. If the points overlap, similarity should be 100%.


title: How to compare the performence between different models?

text: Given a data, I'd like to utilize it by two different models, such as generalized linear model (GLM) and generalized partial linear model (GPLM), and my problem is how to compare the performance of these models? Please give me some ideas! Thanks!


title: in SAS how can I output Number of observation in regression output dataset?

text: I am outputting regression in a dataset but it doesn't contain the number of observations. How can I add this to be added to the dataset? Cheers


title: strucchange breakpoints command: does NA mean no breaks are identified?

text: I am using the package strucchange to analyze a monthly time series. I read it in as a `zoo` object. The series looks something like this:                       date            A            B     2 1979-01-07            0            7               3 1979-01-14            0            4                 4 1979-01-21            0            5                 5 1979-01-28            0            19             6 1979-02-04            0            20      I then run the command:               breakpoints(Conf ~1, breaks=2)      And in some cases, I get the answer:                    Optimal 1-segment partition:           Call:     breakpoints.formula(formula = timeseriesMon_short[, 4] ~ 1, breaks = 2)          Breakpoints at observation number:     NA           Corresponding to breakdates:     NA       Does this simply mean that no structural breaks were identified? I could not find an answer to this in the documentation of the package.


title: Are there some functions in R programming language for generalized partial linear model?

text: I want to use generalized partial linear model to fit a data. So are there some existing functions about GPLM? I'm new about R.


title: Fitting a Poisson GLM mixed model with a random slope and intercept

text: I'm currently working on a series of Poisson time series models trying to estimate the effect of a change in how the counts were obtained (switching from one diagnostic test to another) while controlling for other trends over time (say a general increase in the incidence of disease). I've got data for a number of different sites. While I've been tinkering with GAMs as well, I've fit a series of pretty basic GLMs with time trends in them, then pooling the results. The code for this would look something like this in SAS:               PROC GENMOD data=work.data descending;       model counts = dependent_variable time time*time / link=log dist = poisson;     run;      or this in R:               glm(counts ~ dependent_variable + time + time*time, family="poisson")      Then taking those estimates, and pooling them over the various sites. It's also been suggested to be that I try using a Poisson mixed model with a random slope and intercept for each site, rather than pooling. So essentially you'd have the fixed effect of dependent_variable, then a random effect for the intercept and time (or ideally time and time^2 though I understand that gets a bit hairy). My issue is I have no idea how to fit one of these models, and it seems that mixed models are where everyone's documentation goes suddenly very opaque. Anyone have a simple explanation (or code) on how to fit what I'm looking to fit, and what to look out for?


title: Endogeneity & IV

text: Consider the the following structural model: $y=\beta_1x_1+\beta_2x_2+u$ where $u$ is an iid disturbance term. Suppose $E(u|x_1)=0$ but $E(u|x_2)\neq 0$.For $z_2$ to be a valid instrument, it must satisfy the following two conditions: 1) $E(u|z_2)=0$ 2) $E(x_2|x_1, z_2) \neq E(x_2|x_1)$ **I am struggling to understand the intuition of the relevance condition 2).** Normally for a structural model with one endogenous variable, no $x_1$ just $x_2$, the relevance condition is: 2') $Cov(x_2, z_2)\neq0$ **Why is the condition different when the structural model contains one exogenous variable and one endogenous variable?** Thank-you!


title: How to make a feature vector based on this problem

text: I meet a problem when I wanted to convert raw data into a feature vector in the data collection stage in a data mining project. Let me describe the problem first: Basically, it is a time estimation problem of executing a cluster of blocks. 1.Block has countable number of types such as typeA, typeB, typeC... 2.Each type of block can be more than one in a cluster and a cluster can have different types of blocks. 3.Each block's run time is determined by countable number of factors such as factorA, factorB, factorC... 4.I only have the running time of a cluster of blocks for training and don't have each block's running time. I would like to know how I can make it into a feature vector? The hard part is each block forms something like a structured information and it is hard for me to integrate this structured information into the feature vector.


title: Appropriate negative eigenvalue correction for PCoA of genetic distances

text: I am trying to find the best way to represent genetic distances in a plane so that they may use them as response variables in canonical redundancy analysis (using `rda()` in `vegan`). While there are admittedly many genetic distances to choose from, none of them are Euclidean, which produces negative eigenvalues when we "summarize" them using PCoA (Principle Coordinate Analysis). As far as I know, there are 2 functions to correct negative eigenvalues: Lingoes and Cailliez. The following explanations on how they are computed are pulled out of the `R` documentation for the `pcoa()` function in the `PCNM` package, but they are also well explained in `ape`:   * Lingoes: In the Lingoes (1971) procedure, a constant c1, equal to twice absolute value of the largest negative value of the original principal coordinate analysis, is added to each original squared distance in the distance matrix, except the diagonal values. A new principal coordinate analysis, performed on the modified distances, has at most (n-2) positive eigenvalues, at least 2 null eigenvalues, and no negative eigenvalue.   * Cailliez: In the Cailliez (1983) procedure, a constant c2 is added to the original distances in the distance matrix, except the diagonal values. The calculation of c2 is described in Legendre and Legendre (1998). A new principal coordinate analysis, performed on the modified distances, has at most (n-2) positive eigenvalues, at least 2 null eigenvalues, and no negative eigenvalue. During my explorations, I have not come across these transformations in any population/landscape genetics literature. Many authors use **Fst** -based genetic distances, and favor computing the square root of these distances before applying PCoA, and then removing any remaining negative eigenvalues (these typically account for <1% variance). I have also seen this for **Dps** (proportion of shared allele distance, or Bray-Curtis distance), but in community composition studies. In fact this is a procedure that is recommended in the `CANOCO 5` manual (p.104). While I would normally shrug this off after trying a few a defaulting to what the literature seems to favour, but I explored all of these transformations on a **Fst** -based genetic distance matrix and the downstream results vary considerably. Namely, the genetic variance explained by my models increases by an adjusted R2 of about 0.10 when I use the Cailliez transformation instead of computing `sqrt()` prior to the PCoA. My instinct tells me not to go down that path, but I'm really uncertain why these different corrections yield such drastically different results, and why I should favor the square root correction. Thanks for any insight you may have.


title: Evaluation of probabilistic predictions

text: In the 2010 KDD cup, participants were tasked with estimating the probability that a student would solve a particular exam question. The competition winner was whoever produced the lowest root square mean error (RMSE). Is there some theorem that guarantees that guessing the true probabilities provides the lowest value? Also, it's hard to interpret the results because, the RMSE will be non-zero even if we were perfectly able to guess the true probabilities of a user getting a question right. It would seem that there should be a method of seeing that, say 71% of the questions we marked as the student of having 70% probability of solving were done correctly (and so on for our other predictions), so the algorithm should produce a near zero value


title: Creating consensus from multiple methods of measuring the same entity with some missing values

text: Imagine we have _C_ cars and _D_ drivers, and each driver takes a large subset of these _C_ cars in order to test the rate of fuel consumption for some fixed amount of fuel (let's assume that the number of cars not tested by any given reviewer is relatively small, e.g. <20%). Every driver has their own preferred way of testing, which means that given the same car each driver will get different consumption rates. For example, perhaps Driver 1 takes the cars down a highway as fast as possible, Driver 2 drives through an inner-city route with the headlights on, Driver 3 just sits in the car and turns the air-con to full power, and so on. So we end up with a matrix of fuel consumption values **F** of size _C_ x _D_ with some missing values, and even though the values in each row measure the same underlying notion, there may be a large amount of variation because the testing methods are different. My question is: how can we use this information covering a variety of usage types - but also with missing data - to find a general consensus about which car has best/worst fuel consumption?


title: Is possible use HoG features to train a Neural Network?

text: I was wondering if it is possible use HoG features to train a Neural Network, I know that in the original paper by Dalal and Triggs they used the data generated to train a SVM. If not is possible or _suitable_ do you know a descriptor to get features in order to train a neural network.


title: Measuring the moderating effect in my analysis?

text: I'm trying to analyze if the volume and sentiment (valence) of tweets has an effect on the sales of a product (in my case movies). The DV is sales with IV NumberOfTweets, Valence, DaysSinceRelease. My regression analysis showed that NumberOfTweets is highly significant for Sales (p = .000) while Valence wasn't (p = .731) and, kinda obvious, so was Days Since Release (p = .000). I would now also like to analyze how Valence moderates it effects as my hypothesis states that higher valence leads to higher volume (in turn higher sales). How would ago about measuring this?


title: Does this histogram of residuals indicate that the residuals are effectively random?

text: I am studying a univariate and discrete time series. I know that residuals should be effectively random and have a good fit, and should have a bell shape. **Does the plot below suggest that the residuals are effectively random?** ![Residual Plot](http://hostr.co/file/970/jUNt1ui2JySM/Residui.jpg)


title: A question regarding the central limit theorem

text: I would appreciate if you could please take a look at the following attempt. I am preparing for an upcoming exam, and this is taken from a past exam paper. I would like to know if it is correct. Thank you for your time. Consider the average $Y$ of $n$ independent random variables, each uniformly distributed on $[0,1]$. Find $n$ so that the probability that $Y$ is smaller than $0.51$ is about $0.98$. Make use of the following values of the cumulative distribution function $\Phi$ of the standard normal distribution. $\Phi(0)=0.50000$, $\Phi(0.5)=0.6915$, $\Phi(1)=0.8413$, $\Phi(1.5)=0.9332$, $\Phi(2)=0.9772$. **Attempt** $Y= \dfrac{\sum_{i=1}^{n}Y_{i}}{n}$, where $Y_{i}$, $i=1, 2, ..., n$ is uniformly distributed on $[0,1]$. Then $\mu = E[Y_{i}]=\frac{1}{2}$ and $\sigma^{2}=Var(Y_{i})=\frac{1}{12}$. We have $$Z_{n}=\dfrac{\sum_{i=1}^{n} Y_{i} - n\mu}{\sigma \sqrt{n}}=\dfrac{\sum_{i=1}^{n} Y_{i} - \left(\frac{1}{2}\right)n}{\sqrt{\frac{n}{12}}}.$$ Then $$\dfrac{\sum_{i=1}^{n} Y_{i}}{n}=\dfrac{Z_{n}}{\sqrt{12n}} + \dfrac{1}{2}.$$ Now, by the central limit theorem we have $$P(Y \lt 0.51) = P\left(\dfrac{\sum_{i=1}^{n} Y_{i}}{n} \lt 0.51\right)=P\left(\dfrac{Z_{n}}{\sqrt{12n}} + \dfrac{1}{2} \lt 0.51 \right),$$ rearranging, we get $$P\left(Z_{n} \lt \sqrt{12n} \left(0.51 - \frac{1}{2} \right)\right).$$ we want $P(Y \lt 0.51) \approx 0.98$, then $$P\left(Z_{n} \lt \sqrt{12n} \left(0.51 - \frac{1}{2} \right)\right) \approx 0.98 \Rightarrow P\left(Z_{n} \lt \sqrt{12n} \left(0.51 - \frac{1}{2} \right)\right) \approx \Phi(2),$$ hence $$\sqrt{12n} \left(0.51 - \frac{1}{2}\right) \approx 2.$$ Solving, we get $n \approx 3333.33 $, which is not a whole number, therefore we take $n = 3334$.


title: Find exceptional parameters

text: I have been given an excel document with many rows full of numbers, some rows are marked. Each row represents a case in the clinic, each column represents a research test parameter. I need to find, what's so special about these marks that make them marked, and it has to be in a unified pattern, so I can create an algorithm to find the same exceptions in the Excel file. I'll be using C# to write the program, but I'm sure the question is about what to look into rather than how to actually write the code. I'm brand new to this website, I'm a programmer, not a math genius, so please be patient. I'll check for comments and respond as soon as possible.


title: Automatic keyword generation evaluation

text: I have a simple Text analyzer which generates keywords for a given input text. Till now I have been doing a manual evaluation of it, i.e., manually selecting keywords of a text and comparing them against the ones generated by the analyzer. Is there any way in which I can automate this? I tried googling a lot for some free keyword generators which can help in this evaluation but have not found any till now. I would appreciate any suggestions on how to go about this.


title: Incidence Rate Ratios in R (poisson regression)

text: I would like to run a poisson regression and have R calculate incidence rate ratios, much like the irr option for poisson regression in stata. If anyone can advise how to do this I would be grateful. Thanks PS


title: What is the reason for differences between nbreg and glm with family(nb) in Stata

text: I am a novice Stata user (forced here from R to conform to coauthor's quirky option choices). I need access to the Pearson residuals from a negative binomial regression. I currently have the regression equations specified in this manner:               xi: nbreg cntpd09 logpop08 pcbnkthft07 pccrunion07 urbanpop pov00 pov002 edu4yr ///      black04 hispanic04 respop i.pdpolicy i.maxloan rollover i.region if isser4 != 1,      cluster(state)     est store pd      As I understand it, to get the Pearson residuals I need to undertake the same regression using the glm function. Something like this:               xi: glm cntpd09 logpop08 pcbnkthft07 pccrunion07 urbanpop pov00 pov002 edu4yr ///      black04 hispanic04 respop i.pdpolicy i.maxloan rollover i.region if isser4 !=1, ///      family(nb) vce(cluster state)      Some details to note: 1)I use an if statement to select a subset of the total records, 2)I am using robust standard errors clustered by "state" Sorry for not including a replicable example at this stage, but I am new enough to Stata that I am not sure what each part of these statements is doing exactly and wanted the formatting to be exact. Is there something basic in the coding that would explain differences in coefficients and standard errors between these methods as I have them written out? Should I be worried about these differences? How can I make the glm version conform to the nbreg version? A secondary question since I have your attention already, how do I save, aggregate, and export the Pearson residuals assuming I get the glm function to work properly and am able to run it on each of the 6 regression models in our paper (Each model has the same observations, just different dependent and independent variables)? Ideally I would have a table with six columns of residuals in any format readable by R (including dta or other Stata defaults). Thanks in advance, Chris


title: Kenward-Roger correction applied to $t$-test in LME

text: Kenward and Roger (1997) proposed an approximation for an $F$-test of fixed effects in a linear mixed-effects model through (1) inflating the estimated variance-covariance matrix of the fixed and random effects and (2) adjusting the degrees of freedom (DF). In their paper they state that, when the approximation is applied to a $t$-test (a special case of the $F$-test), it is "comparatively straightforward and essentially recovers Satterthwaite's (1941) approximation". My understanding is that Satterthwaite's method involves only adjustment of degrees of freedom. My question is, when applying the Kenward- Roger correction to a $t$-test for a fixed effect, would we only need DF adjustment, or would we need correction for both DF and the variance- covariance structure? Thanks!


title: Why use a certain measure of forecast error (e.g. MAD) as opposed to another (e.g. MSE)?

text: MAD = Mean Absolute Deviation MSE = Mean Squared Error I've seen suggestions from various places that MSE is used despite some undesirable qualities (e.g. http://www.stat.nus.edu.sg/~staxyc/T12.pdf, which states on p8 "It is commonly believed that MAD is a better criterion than MSE. However, mathematically MSE is more convenient than MAD.") Is there more to it than that? Is there a paper that thoroughly analyzes the situations in which various methods of measuring forecast error are more/less appropriate? My google searches haven't revealed anything. A similar question to this was asked at http://stackoverflow.com/questions/13391376/how-to-decide-the-forecasting- method-from-the-me-mad-mse-sde, and the user was asked to post on stats.stackexchange.com, but I don't think they ever did.


title: Is the p-value still uniformly distributed when the null hypothesis is composite?

text: When the null hypothesis is simple (i.e., has only distribution of the sample) and the only distribution of the sample is continuous, the p-value can be shown to be uniformly distributed over (0,1). When the null hypothesis is composite (i.e., has more than one sample distributions), is the p-value still uniformly distributed over (0,1)? Why?


title: Nonlinear regression: Confidence intervals on transformed or untransformed parameters?

text: Suppose I am using a standard inhibition model to find biochemical parameters that fit my data. The equation is: $y = \frac{A}{{1 + \exp \left( {\ln \left[ S \right] - \ln IC_{50}} \right)}} $ where $\left[ S \right]$ is the concentration of my inhibitor and $IC_{50}$ is the concentration of my inhibitor at which the measurement (with a maximum $A$) is reduced by half. Which of these approaches should I take?   1. enter the equation into the NLR procedure as-is (FYI: I am using SPSS) and let it fit the values of $A$ and $IC_{50}$ in the same manner, i.e. with all of the assumptions of OLS regression parameters. Or...   2. enter the equation as  $y = \frac{A}{{1 + \exp \left( {\ln \left[ S \right]- \ IC_{50}^* } \right)}} $ where $ IC_{50}^* = \ln IC_{50}$. This of course would require a transform on the output parameter and confidence limits of $ IC_{50}^*$ giving me asymmetric error bars. Which of these strategies is most rigorous? My instincts suggest the 2nd option, as $ IC_{50}$ is actually bound by 0 and is thus more likely to be log-normal rather than normal. Any help (direct answers, references, etc) is appreciated.


title: Is this an example of Pooled OLS on Panel Data?

text: I am looking at a study that analyzes the effect of an infrastructure index on infant mortality and child mortality rates. The database has (asset) quintile level data for 47 different countries (from surveys that used the same methodology in all countries, althought they took place @ different years), thus there are 47*5=235 observations in total. The authors control for some country level factors (i.e. do not vary among quintiles) such as as income per capita and national literacy rate of females. All other independent variables are country quintile specific. They then give these results. They say they calculated the regressions using Seemingly Unrelated Regressions (SUR), but they just look like usual OLS regressions to me. I have the following questions: 1) By using SUR, are they simply running a pooled ols over all of the 235 data points? The authors say they are using "cross country quintile level data", and simply state that they used SUR for regressions. 2) I assume the approach 1) would be incorrect since a commentator for the paper suggested setting up the data as a Panel formed over quintiles and countries; and then using FE instead to account for country level fixed effects. I've never heard of a Panel formed this way, does this make sense? And would the FE method makes sense then? If anyone's interested the study in question is: http://elibrary.worldbank.org/doi/pdf/10.1596/1813-9450-3163


title: How to calculate the variance for mean of means?

text: I have been reading up on how to calculate variance for different situations and the following has me confused….. If I take a number of samples (e.g. choose 5 kids in a class and ask them how many pets they have), each time I go into a class and choose 5 kids I am going to get a different mean number of pets, and then each sample will have its own variance. If I then took the mean of all my samples and added them up and then divided them by the number of the sample I had I'd get the mean of the means. The part I don't understand is how do I calculate the variance of the mean of means ? Do I subtract the each mean from the mean of means, square the answer and then divide by n ?


title: Sampling in $\mathbb{R}^n$ with roughly equal Voronoi cells

text: N random points in a ball in $\mathbb{R}^n$ induce a Voronoi splitting or tessellation.   Is there a way of random-generating points so that the volumes of the Voronoi cells are roughly equal ? This is 3 related questions, for $L_1, L_2$ and $L_\inf$ norms (in 3d, random splits of an octahedron, sphere, cube). (Where does this come from ? In searching for a minimum of a noisy function of $n$ variables, one may start N parallel searches from N random start points, e.g. 10 start points in 5d.   Small Voronoi cells are then inefficient, equal-sized cells most efficient.)


title: Is there a method to estimate distribution parameters given only quantiles?

text: is there a way to fit a specified distribution if you are only given a few quantiles? For example, if I told you I have a gamma distributed data set, and the **empirical** 20%, 30%, 50% and 90%-quantiles are, respectively:                     20%       30%       50%       90%      0.3936833 0.4890963 0.6751703 1.3404074       How would I go and estimate the parameters? Are there multiple ways to do that, or is there already a specific procedure? **more edit:** I don't specifically ask for the gamma distribution, this was just an example because I worry I can't explain my question appropriately. My task is that I have some (2-4) given quantiles, and want to estimate the (1-3) parameters of a few distributions as "close" as possible. Sometimes there's an (or infinite) exact solution(s), sometimes not, right?


title: Is positive coefficient of price correct in a multiple regression model

text: I am currently undertaking forecasting of energy sales (kWh) for our industrial customers. From historical data gathered from 1993 to 2013, a graph of price per kwh against sales kwh shows a positive relationship. Going ahead to generate the model yielded a positive coefficient of price with energy sales. The other independent variables are GDP of the industrial sector and the lag of sales for 1, 4 and 6. What can explain the positive coefficient of price with sales? Plotting GDP of the industrial sector yields a positive relationship with sales (kWh). Can someone help? I need to understand why we have a positive slope of price with sales. The historical data clearly confirms that.


title: Can I delete excessive number of multivariate outliers, like over 10% in sample?

text: I have a dataset with around 9000 cases, I'm running a factor analysis and I have found that 1100 cases are identified as being a multivariate outlier. Is it alright for me to got ahead and delete it?


title: How to customize a link function to perform a logistic regression?

text: My data was collected using Randomized Response Technique. So I have additional variability into the data. I have a binary response variable. Should I customize a logit link function to incorporate the known probabilities of Randomized Response Technique to perform a logistic regression? I would like to know if this approach is appropriate or how one might correctly perform a regression with these kinds of responses. If a customized link function is appropriate, I should customize the logit $\log \left({p \over 1-p}\right)$ to the customized $\log \left({p-0.25 \over 0.75-p}\right)$. Can anyone help me to do that? I don't have too much experience using R.


title: Regression with covariates on time-series data

text: I have monthly data at physician level (sales and marketing activity data for last 12 months for 3000 physicians). I want to find out the impact of marketing activity on sales. There are other covariates and independent variable is continuous. How should I approach the model?


title: scale of measurements in sign test and signed rank tesr

text: Can some one help me to understand the right scale of measurement to use in sign test as well as in signed rank test? I read a certain book it says nominal scale is used in sign test while ordinal scale is used in signed rank test, the other online material says ordinal scale is used in both sign test and signed rank test im confused.


title: Record Linkage Using Fellegi-Sunter Model

text: I am trying to create a record linkage system using the fellegi-sunter model.I am following this paper http://digital.library.okstate.edu/etd/SHIN_okstate_0664M_10668.pdf. I am not understanding clearly how to calculate the upper and lower threshold after calculating the M,U probabilities. The paper states that the threshold values depend on false positive and false negative values.But how to calculate the false positive and false negative values?


title: A question about warning that performing k-fold CV with caret

text: I am trying to `createFolds` function in `caret` to use k-fold cross- validation in R. But I came across this warning:                Warning messages:     1: In foldVector[which(y == dimnames(numInClass)$y[i])] <- sample(seqVector) :       number of items to replace is not a multiple of replacement length      Has anyone come across this before & solved it? I am using this R code:               > flds<-createFolds(dataset,k=5, list=TRUE, returnTrain=FALSE)      The data-set includes my response variables that I am going to use with Linear Discriminant Analysis > str(dataset) 'data.frame': 200 obs. of 162 variables: $ V2 : Factor w/ 201 levels "6.74361083760967",..: 150 190 108 157 180 7 139 116 175 49 ... $ V3 : Factor w/ 201 levels "1.88828398539874",..: 114 177 42 155 143 173 144 81 142 117 ... [list output truncated]


title: Under-forecasting in Regression

text: I have to do forecasting of sales that is how much sales of a product is going to happen in a particular store. I have time series data for last two years and doing forecasting for 2014. The variables are promotion flag ( Yes/ No ), promotion period, location in a store, price discount. These all are categorical variables. For this I am using regression method where, dependent variable is sales, and independent variables are categorical variables mentioned above. This analysis is done in SPSS where I have used step-wise and backward regression. I want to know, the regression model is under-forecasting? Is there a way to improve the forecast?


title: How to include noise in clustering evaluation?

text: When evaluating clustering methods which do have a definition for noise points (like dbscan), how noise will affect evaluation? Consider a clean dataset like well known Iris dataset. There is no noise in it. In evaluation metrics like homogeneity should we remove noise labels and then compute the metric? Ex:               True_labels = [1, 1, 1, 2, 2, 1]     Expected_labels = [a, a, -1, b, -1]      Do we need to alter these labels some how like this?               True_labels = [1, 1, 2]     Expected_labels = [a, a, b]


title: prediction interval of a new prediction out of a weighted linear regression model

text: This question is based on the example 11.1 out of the book Applied linear statistical models of Kutner, Nachtsheim, Neter and Li. You can find the data here. First they calculate a simple linear model:               Call:     lm(formula = Bloodpressure ~ Age, data = data)          Residuals:          Min       1Q   Median       3Q      Max      -16.4786  -5.7877  -0.0784   5.6117  19.7813           Coefficients:                 Estimate Std. Error t value Pr(>|t|)         (Intercept) 56.15693    3.99367  14.061  < 2e-16 ***     Age          0.58003    0.09695   5.983 2.05e-07 ***     ---     Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1          Residual standard error: 8.146 on 52 degrees of freedom     Multiple R-squared:  0.4077,    Adjusted R-squared:  0.3963      F-statistic: 35.79 on 1 and 52 DF,  p-value: 2.05e-07      But the residual plots show that the variance of the error term increases with predictor `Age`. Long story short, they use weights to account for this heteroscedasticity and they get:               Call:     lm(formula = Bloodpressure ~ Age, data = data, weights = weights)          Weighted Residuals:         Min      1Q  Median      3Q     Max      -2.0230 -0.9939 -0.0327  0.9250  2.2008           Coefficients:                 Estimate Std. Error t value Pr(>|t|)         (Intercept) 55.56577    2.52092  22.042  < 2e-16 ***     Age          0.59634    0.07924   7.526 7.19e-10 ***     ---     Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1          Residual standard error: 1.213 on 52 degrees of freedom     Multiple R-squared:  0.5214,    Adjusted R-squared:  0.5122      F-statistic: 56.64 on 1 and 52 DF,  p-value: 7.187e-10      Now my question is, how do I calculate the prediction error if I get a new patient with `age=25`? Is it correct that the prediction interval will be smaller for lower ages?