title: What are the benefits of normality assumption for AR(I)MA models?

text: I know normality assumption is not necessary for all of the ARIMA models. My question is that if we have a non-normal time series, is it better to transform it to normal state by transformations like Box-Cox and then use it in the model? what are the pros and cons of these transformation in ARIMA models?


title: SPSS and Stata output different

text: I'm Stata-proficient but learning SPSS for my new position. I am using a simple dataset to do very basic regressions and comparing to see if the results are the same. They're not. I'm close, but the magnitudes of the betas and significance are slightly different. The data was copy and pasted into each from an Excel; I didn't use a Stata file in SPSS, or vice versa. For SPSS, I did not weight it, it's using listwise deletion, and it's on the "enter" method. I presume Stata is doing the same, as its default (but correct me if I'm wrong and it's a different default!). Any ideas on what else to check? I'm doing just a simple linear regression. Data: https://www.dropbox.com/s/8g31cjf8vr69i44/rwj%20county%20data.xls?dl=0 ### Syntax SPSS               REGRESSION     /MISSING LISTWISE     /STATISTICS COEFF     /DEPENDENT FreeLunch     /METHOD=ENTER FoodInsecure Rural Female @18 Hispanic.      (Or, for point and click, analyze->regression->linear; forces a choice under "method" for enter/stepwise/remove/backward/forward.) Stata               reg percentfreelunch percentfoodinsecure rural female under18 hispanic      Data is in Excel and was pasted into both. ### Results SPSS                          Var. |   Unst.B | Std.Err. |   St.B |     t | Sig.          (Constant) | -139.616 |   66.652 | -2.095 | .045     % Food Insecure |    2.785 |     .674 |   .546 | 4.131 | .000               Rural |     .131 |     .048 |   .404 | 2.701 | .011              Female |    2.657 |    1.170 |   .372 | 2.270 | .031                < 18 |    -.416 |     .583 |  -.145 | -.715 | .480            Hispanic |    1.156 |     .236 |  1.092 | 4.905 | .000      Stata                              Var. |     Coef. | Std.Err. |     t | P>|t|     percentfoodinsecure |   2.76532 | .6741544 |  4.10 | 0.000                   rural |  .1378976 | .0495354 |  2.78 | 0.009                  female |  2.826711 | 1.204272 |  2.35 | 0.026                 under18 | -.3799895 |  .588423 | -0.65 | 0.523                hispanic |  1.168375 | .2398765 |  4.87 | 0.000                   _cons | -149.3858 |  69.0891 | -2.16 | 0.039


title: Fitting Model with MCMC?

text: In fitting a model with several variable I have found extremely useful a method involving the minimization of the Chi Square using a MCMC approach. In particular, I followed this tutorial http://sciencehouse.wordpress.com/2010/06/23/mcmc-and-fitting-models-to-data/ However, I cannot find around any other source. Since I will use this method for a scientific work, I would like to a) be sure that the method is theoretically correct b) cite a peer-reviewed source. Running several simulation with complex model (6/7 parameters) showed me that the model actually works really fine, retrieving the parameter within a small error most of the time. In particular, I am not sure about the use of $exp(-\chi^2) $ for the likelihood. Is this justified? I also tried to calculate the ratio by just dividing the oldChi2 with the newChi2, and the results were, again, satisfying. Is there any relevant difference between these two methods? Thank you


title: Understanding the output of lme

text: I am having some difficulties interpreting the results of an analysis perfomed using lme. I conducted an experiment where the subjects had to estimate the time elapsed in a task involving a spatial measure (e.g. subjects watched a video game where a car travels a certain distance). My goal is to determine if there is a linear relation between perceived time and the space traveled by the car. Subject underwent 2 different conditions (e.g. the car was traveling two distances), and each condition was repeated twice. Therefore, given the nature of the experimental design involving repeated measures, I cannot use a simple linear model but I have to use a mixed effects linear model with a random intercept for subject. I use the R language, and I adopted this formula to solve my problem               library(nlme)     summary(lme(Time ~ Distance, data = my_Table, random = ~1 | Subject))      The output that I get is:               Linear mixed-effects model fit by REML      Data: my_Table            AIC      BIC    logLik       608.315 618.4454 -300.1575          Random effects:      Formula: ~1 | Subject             (Intercept) Residual     StdDev:    2.964139 4.919044          Fixed effects: Time ~ Distance                                               Value Std.Error DF  t-value p-value     (Intercept)                            5.518714 0.8212930 64 6.719543  0.0000      Distance                              0.013092 0.0053225 64 2.459718  0.0166      Correlation:                                             (Intr)      Distance                              -0.415          Standardized Within-Group Residuals:            Min         Q1        Med         Q3        Max      -1.2904862 -0.4302117 -0.2593047  0.2081889  5.0914269           Number of Observations: 95     Number of Groups: 30       Now, what is this output telling me? As far as I understand from studying the literature that I found online (I am not a statistician...), if I am rght these results are telling me that there is a linear relation between time and space, since b = 0.013092 and p-value is significant. Now, if the p-value was not significant, this would have meant that there is no linear relation between time and space?


title: Calculating necessary sample size using bootstrap

text: I have recently come across a website (http://www.surveysystem.com/sscalc.htm) that returns the sample size given the following inputs: confidence level, confidence interval, and population. I assume this is done by rearranging a hypothesis test under a CDF, I guess using the standard normal distribution? But, if one does not believe the data reflects this particular CDF, how would you go about using bootstrapping to arrive at a version of the sample size that is data driven and not limited to the ~N(0,1) distribution? Further, I am interested in this procedure for stratifications of a sample. Help on either, especially the first paragraph, would be appreciated.


title: ROC curve for discrete classifiers like SVM: Why do we still call it a "curve"?, Isn't it just a "point"?

text: In the discussion : how to generate a roc curve for binary classification, I think that the confusion was that a "binary classifier" (which is any classifier that separates 2 classes) was for Yang what is called a "discrete classifier" (which produces discrete outputs 0/1 like an SVM) and not continuous outputs like ANN or Bayes classifiers ... etc. So, the discussion was about how the ROC is plotted for "binary continuous classifiers", and the answer is that the outputs are sorted by their scores since the outputs are continuous, and a threshold is used to produce each point on the ROC curve. My question is for "binary discrete classifiers", such as SVM, the output values are 0 or 1. So the ROC produces just one point and not a curve. I'm confused as to why we still call it a curve?!! Can we still talk about thresholds? How can one use thresholds in SVM in particular? How can one compute the AUC?, Does cross-validation play any role here? thanks in advance.


title: Why is random assignment important in stratified sampling?

text: **Background** I raised this question because of an argument I am having with a question from user697473 here. The title of his question is "Formal definiton of random assignment." In the post he makes a claim that one can get unbiased estimates of the difference of treatment effect even if one group that is represented in the population is not sampled. At least that is the way I interpret the claim. He points to a paper coauthored by Don Rubin on causal modeling here. Both he and StasK argue through examples that this claim is true. This issue was also raised in another question that he posted here. In that post titled "Random assignment: why bother?" gung expresses skeptism. In the first question my answer included a statement to the effect that without sampling all possible assignments of a potential confounding factor you cannot construct an estimate of the difference in treatment effect without assigning samples to each stratum for the confounding factor. **The Issue** I am confused about the claim. As I interpret the claim I am sure it is false. The claim is supported by two examples one given by user697473 and the other by StasK. The examples are confusing to me and in neither case is a proof or demonstration given to show that the claim is true. I could be misinterpreting the claim but based on my interpetation it is false. To illustrate my point let's look at a famous example from the design of experiments book of the Godfather of randomization Sir Ronald Aylmer Fisher: "The Lady Tasting Tea". In case you don't know the example I will describe it and provide a modification to clarify the example. The lady claims that she has the ability to determine just by tasting a cup of tea whether the tea was poured first or the milk was poured first. Fisher poses an experiment to test her claim by providing some cups with tea poured first and others with the milk poured first. If the Lady has the expertise that she claims she should be able to correctly identify the types better than a person who guesses at random. So the experiment is design to see if her probability of correct classification is better than 0.5. So he randomly provides the Lady with cups from both groups (in the example there were 6 cups 3 with milk poured first and 3 with tea poured first. Now let me modify the experiment slightly. Suppose I have three brands of tea and for my population these brands are served at a tea parties equally. The question I want to ask is when given a cup of tea at a tea party can the Lady predict whether milk was poured first or the tea was poured first. I want to test whether her prediction accuracy is greater than 0.5 in the setting of tea parties. So I will apply Fisher's experiment to estimate this probability when the Lady is at one of these tea parties.   **Note:** I do not care whether or not she can differentiate the brands. Suppose that in the population of tea parties the Lady's ability is: **Brand A** Can predict tea first with probability 0.7 and can predict milk first with probability 0.7. **Brand B** Can predict tea first with probability 0.75 and can predict milk first with probability 0.75. **Brand C** Can predict tea first with probability 0.5 and can predict milk first with probability 0.6. So she does no better than chance with Brand C but can do better than chance with Brands A and B. For the experiment I only provide Brands A and B for her to taste. All I want to know is her ability to predict correctly in the normal tea party situation. **My Questions**   1. Can I get an unbiased estimates of her prediction capabilities with this experiment?   2. If the answer to (1) is no, is user697473's claim false or have I misinterpreted it? Randomization is used to avoid bias that is produced by confounding variables. This is important when trying to draw inference from a sample to a population. In my example the Brands are confounders. The treatment is tasting the cups. In the population she would get Brand A, Brand B and Brand C each 1/3rd of the time. If I know this and I sample Brands with unequal probabilities I claim that I can get an unbiased estimate of her prediction capabilities by taking a specific weighted average of the estimates of the proportions if I use a stratified random sample with strata totals n$_1$, n$_2$ and n$_3$ all greater than 0. But I cannot if n$_3$ = 0. Furthermore, there is not other estimate I can calculate from a sample where allocation to Brand C equal 0 that will be an unbiased estimate.


title: How to interpret the odds ratio in a logistic regression with proportion as a response variable

text: I have a glm model for some data with a proportion as the outcome variable as follows:               xi: glm pos_tests i.covariate_1 covariate_2 , family(binomial total_tests)     link(logit) vce(robust) eform      In this model, `pos_tests` is the number of positive laboratory results and `total_tests` in the number of all tests conducted. `i.covariate_1` and `covariate_2` are categorical and continuous independent variables respectively.                     pos_tests  | Odds Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]         -------------+----------------------------------------------------------------        _Icovariate_1 |   2.535989   .9193459     2.57   0.010     1.246149    5.160891          covariate_2 |    .976045   .0270064    -0.88   0.381     .9245231    1.030438      My question is how do I interpret the significant result for covariate_1. Does it mean that the odds of all test results being positive are 2.54 time greater for those individual with covariate_1 or is there some other way of interpreting. Thanks.


title: Is Matlab/octave or R better suited for monte carlo simulation?

text: I started to do Monte Carlo in R as a hobby, but eventually a financial analyst advised to migrate to Matlab. I'm an experienced software developer. but a Monte Carlo beginner. I want to construct static models with sensitivity analysis, later dynamic models. Need good libraries/ algorithms that guide me. To me seems that R has excellent libraries, and I suspect mathlab is preferred by inexperienced programmers because of the easy pascal-like language. The R language is based on scheme and this is hard for beginners, but not for me. If Matlab/ Octave has not advantages on the numerical/ library side I would stick with R.


title: How to report a permutation test?

text: I want to compare a parameter between two independent groups of volunteers (see MWE below). All volunteers had baseline reading (V1) and a subsequent reading by a later visit (V2) for intervention (`Interv`; Yes/No). The parameters are gene expression data and expressed in fold change (i.e. V2/V1). For brevity, I mentioned only one `Par1`. I came from here, but I didn't find what I want though how to report was updated in the question but it was not addressed yet I guess, so it might deserve a separate question like this one here. To report the central tendency, I am intending to use the `median` as it is less affected by outliers, may be with `MAD` (Median Absolute Deviation), though not very often used, but I think is the equivalent to `SD` for the `mean`. Subjects without intervention were (`n=15`), whereas those with intervention were (`n=19`). I am not a statistician, but I know this from my field that when it comes to gene expression data it is always good to log transform the data before doing analysis (to make the data mimic the normal distribution). So to report the `median` I used the fold change reading because it is more intuitive to the reader and makes more sense than reporting log values. To measure the difference between the two groups as inferential statistics, I am more inclined to use the `permutation` techniques available in `R`, like the `coin` package, than to use `t-test`, thinking that is more robust and less demanding of assumptions (confession: biologists do not respect these aliens, they will meet anybody except assumptions). After log transformation, I did permutation test, wilcox test (Mann-Whitney-U- test), and t-test just for sake of comparison and learning.   The dataset (use `dget(<code below>)` to get the dataset:               > dput(dfFold)     structure(list(Interv = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L,      1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,      2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("Yes",      "No"), class = "factor"), variable = structure(c(1L, 1L, 1L,      1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,      1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("Par1",      "Par2"), class = "factor"), Fold = c(181.552149430047, 365.287611324048,      1.37615543475921, 29.4749766533613, 7.30460900092265, 0.55729573732714,      0.856344617183982, 12.2151438678328, 8.22023432258123, 7.52091826273063,      11.3960901635291, 141.40637778502, 1.19646833351009, 4234.12741125762,      9410.18608208175, 0.11188317488791, 3.39142462873992, 1.54986235136038,      0.138010480189387, 0.150720706883034, 8.58182709551405, 1.65995914129576,      1.02664060360863, 40.9872498215584, 51.7412085216083, 3.64546611583102,      5.64546161492277, 8.92993632885674, 4.71129389070826, 1497.44850394305,      0.0209326180515642, 2.90945271783875, 5.09753127386217, 9.85421224939179     ), Log = c(5.20154293711853, 5.90068501952248, 0.31929369443907,      3.38354165434038, 1.98850551891357, -0.584659233229834, -0.15508239368987,      2.50267648260455, 2.10659871505853, 2.01768823991267, 2.43327032846195,      4.95163785700696, 0.179374162082163, 8.35093254365055, 9.14954800731289,      -2.19030003343294, 1.22125007758044, 0.438166121413757, -1.98042565330713,      -1.89232677812372, 2.14964683902473, 0.506792988388275, 0.0262919218996518,      3.71326103837802, 3.9462545340732, 1.29348423542785, 1.73085196838103,      2.18940926483575, 1.54996258170508, 7.31151794138033, -3.86644666427115,      1.0679649940249, 1.62875635858491, 2.28789900330443)), .Names = c("Interv",      "variable", "Fold", "Log"), class = "data.frame", row.names = c(NA,      -34L))      How it looks like:               > head(dfFold)       Interv variable        Fold        Log     1    Yes     Par1 181.5521494  5.2015429     2    Yes     Par1 365.2876113  5.9006850     3    Yes     Par1   1.3761554  0.3192937     4    Yes     Par1  29.4749767  3.3835417     5    Yes     Par1   7.3046090  1.9885055     6    Yes     Par1   0.5572957 -0.5846592      To report the central tendency (I included `mean` also):               > median(dfFold[dfFold$Interv=="Yes",]$Fold,)     [1] 11.39609     > median(dfFold[dfFold$Interv=="No",]$Fold,)     [1] 3.645466     > mean(dfFold[dfFold$Interv=="Yes",]$Fold,)     [1] 960.8452     > mean(dfFold[dfFold$Interv=="No",]$Fold,)     [1] 86.71587      Now the the tests, `permutation`, then `wilcox` test from `coin`, and `t-test` from `stats` package:               > pvalue(oneway_test(Log~Interv, data=dfFold ,distribution=approximate(B=9999),conf.int=TRUE))     [1] 0.03590359     99 percent confidence interval:      0.03128298 0.04097196           Warning message:     In independence_test.IndependenceProblem(object, check = check,  :       additional arguments conf.int will be ignored     > set.seed(1234)     > pvalue(wilcox_test(Log~Interv, data=dfFold,distribution=approximate(B=9999),conf.int=TRUE))     [1] 0.05890589     99 percent confidence interval:      0.05300214 0.06523398           > t.test (Log~Interv, data=dfFold)              Welch Two Sample t-test          data:  Log by Interv     t = 2.1519, df = 27.43, p-value = 0.04036     alternative hypothesis: true difference in means is not equal to 0     95 percent confidence interval:      0.09772892 4.04392270     sample estimates:     mean in group Yes  mean in group No               3.183037          1.112211      I got a warning message from the `oneway_test()`, anybody can help explain that?   So permutation test with approximation (B=9999) the difference was significant at 0.05 level, in agreement with the `t-test`, but not with the `wilcox` test. **Now the question:**   If I were to report this difference based on the permutation test what should I mention? Do I have to include the confidence intervals of 99% of `oneway_test()`? (see the output), how to change it to 95%? and what is the difference in concept of statistic between `oneway_test()` and `t-test`? isn't `permutation` a computer-intensive technique that can be used with _many_ statistical tests? but when I do it with `oneway_test()` what kind of test this function eventually uses? Any help to raise my statistical knowledge would be appreciated. Also few lines of comment on the result as one would be expected to write in a report would be much appreciated alike.


title: Estimation of parameters via two component models

text: Given a model with two coponent expression: $$G_R(w) = \frac{k^2 w^{2\alpha} + k w^{\alpha}cos(\frac{\alpha \pi}{2}) }{1+k^2 w^{2\alpha}+2 k w^{\alpha}cos(\frac{\alpha \pi}{2})}$$ $$G_I(w) = \frac{ k w^{\alpha}sin(\frac{\alpha \pi}{2}) }{1+k^2 w^{2\alpha}+2 k w^{\alpha}cos(\frac{\alpha \pi}{2})}$$ where $p = (k, \alpha)$,are model parameters which are needed to fit via experimental data $(w, G^*(w))$. The complex $G^*(w)$ are expressed as $G^*(w) = G_Rdata(w) + iG_Idata(w)$. One standard method for estimation of parameters is nonlinear least square method by minimizing error $$E(p)=\sum_{i=1}^{N}\left| G_R({w}_{i};p) -{ G_Rdata}_{i}\right|^2 + \left| G_I({w}_{i};p) -{ G_Idata}_{i}\right|^2 $$ .However, there is another routine to estimate model parameters, say, we only need one model equation,e.g., by minimizing error $$\tilde{E}(p)=\sum_{i=1}^{N}\left| G_R({w}_{i};p) -{ G_Rdata}_{i}\right|^2 $$ to get the estimation of parameters, and most of the case, we produce an excellent fit to the real component part of data, and provide a complimentary prediction of $G_I({w}_{i};p)$. If we have a poor fitting via second method, then it will not improve with the first one. Could we say the two approaches for estimation of parameters have equivalent performances? Is there any specific mathematical theory on this issue?


title: Upper and lower limits for numbers of runs

text: On this website about statistical process control there is a chart which gives the upper and lower limits for the number of runs in a set of data for the data set to be considered taken from a process that is showing only "common cause" variation. My question: is there any theoretical justification for the limits given, or are the limits context dependent? If there is a justification I would appreciate links to further reading.


title: Econometrics: multivariate cointegration techniques

text: Dependent    Y         R/E        N         W       D98       CONSTANT       G          1.143     -O.656     1.029      0.173    0.104      17.31                   (2.757)   (-2.432)   (1.896)    (0.896)  (2.126)    (3.335)           Adjusted R-squared=0.59      F-statistic=3.218[0.045]      Jarque-bera statistic(with 8 degrees of freedom)=18.05[0.021]      Notes: G=government spending, Y=national output, R/E=tax share(revenue/expenditure), N=population, W=wage rate and D98=1998 dummy; t-statistics in parentheses and p-value in brackets Comment on the overall significance of the model and interpret the estimated coefficients


title: How to apply LASSO to IRLS (logistic regression)?

text: I have programmed a logistic regression using the IRLS algorithm. I would like to apply a LASSO penalization in order to automatically select the right features. At each iteration, the following is solved: $$\mathbf{\left(X^TWX\right) \delta\hat\beta=X^T\left(y-p\right)}$$ Let $\lambda$ be a non-negative real number. I am not penalizing the intercept as suggested in _The Elements of. Statistical Learning_. Ditto for the already zero coefficients. Otherwise, I subtract a term from the right-hand side: $$\mathbf{X^T\left(y-p\right)-\lambda\times \mathrm{sign}\left(\hat\beta\right)}$$ However, I am unsure about the modification of the IRLS algorithm. Is it the right way to do? * * * Edit: Although I was not confident about it, here is one of the solutions I finally came up with. What is interesting is this solution corresponds to what I now understand about LASSO. There are indeed **two** steps at each iteration instead of merely one:   * the first step is the same as before : we make an iteration of the algorithm (as if $\lambda=0$ in the formula for the gradient above),   * the second step is the new one: we apply a soft-thresholding to each component (except for the component $\beta_0$, which corresponds to the intercept) of the vector $\beta$ obtained at the first step. This is called Iterative Soft-Thresholding Algorithm. $$\forall i \geq 1, \beta_{i}\leftarrow\mathrm{sign}\left(\beta_{i}\right)\times\max\left(0,\,\left|\beta_{i}\right|-\lambda\right)$$


title: How to determine the x where P(X > x)=.67 from a normal distribution?

text: I wish to calculate the value on a normal distribution (having the mean and stdev), where P(Z >= ?) = 67%. I'm new to stats, so in case my question is unclear, I will draw a picture: ![enter image description here](http://i.stack.imgur.com/xQqCm.png)


title: Binomial test for a binary variable

text: I have a binary variable (which takes values 0,1). I have about 100k records of it. How do I determine if it follows the binomial distribution? (I'm bascially trying to test for normality. And, if the data is not normal, I might have to apply a transformation to get the variable into a binomial distribution.) * * * Hey, thank you folks for clearing this up. This was an effort as a prelude to cluster analysis. I also understand that normality of variables is more a nice to have feature for cluster analysis and that the distance measures would be valid even otherwise. Your views?


title: How do I perform a regression on non-normal data which remain non-normal when transformed?

text: I've got some data (158 cases) which was derived from a Likert scale answer to 21 questionnaire items. I really want/need to perform a regression analysis to see which items on the questionnaire predict the response to an overall item (satisfaction). The responses are not normally distributed (according to K-S tests) and I've transformed it in every way I can think of (inverse, log, log10, sqrt, squared) and it stubbornly refuses to be normally distributed. The residual plot looks all over the place so I believe it really isn't legitimate to do a linear regression and pretend it's behaving normally (it's also not a Poisson distribution). I think this is because the answers are very closely clustered (mean is 3.91, 95% CI 3.88 to 3.95). So, I am thinking I either need a new way of transforming my data or need some sort of non-parametric regression but I don't know of any that I can do in SPSS.


title: What is the normality test for binary data?

text: While most of the tests for normality help for continous variables, is there a way to test the normality assumptions for binary variables. From what i've read on wiki, the K-S test can be applied for continous variables. How do we perform normality tests for binary (or even categorical variables)?


title: Factors missing in summary and Anova output from glm in R

text: I'm running a GLM model in R, but when I add additional variables, the degrees of freedom (df) for two of the variables (level and type, see models below) goes from 1 to 0, and because of this there is no information about them in either the Anova table or the summary output. Interactions involving these two variables however, continue to have one degree of freedom. I can't figure out why the DF goes to zero when I add additional covariates, or what I can do to address this and was hoping someone could point me in the right direction. The code is below. To try to figure out what is going on, I have run three models:   1. Model 1 with just a few factors (plus their interactions)   2. Model 2 with all of the factors (but no interactions)   3. Model 3 with all of the factors plus the interactions I want The first model runs fine, with the two variables level and type having one df (and therefore output exists for these two variables in the Anova and summary outputs). But for the second and third models, the two variables level and type each have zero df (although each of their interactions with medium has one df).               modeloutput <- glm(model, data=dataset, family=binomial)     ## formulas:     # Model 1:     model <- as.formula(success ~ 1 + medium + level + type + motivation + level*medium +                        type*medium + motivation*medium)     # Model 2:     model <- as.formula(success ~ 1 + course + medium + level + type + motivation +                          ethnicity + gender + age + enrollment + Pell + TANF + GPA +                          onlineExp)     # Model 3:     model <- as.formula(success ~ 1 + course + medium + level + type + motivation +                          ethnicity + gender + age + enrollment + Pell + TANF + GPA +                          onlineExp + level*medium + type*medium + motivation*medium)      I feel like I must be missing something obvious, but I'm just not seeing it. If anyone has any tips about what I can do to find out more information about why this is happening (and about what my options are in trying to address it), I would be very grateful for the help! Edit to include sample Anova output:                                  LR Chisq Df  Pr(>Chisq)     course             339.1891235  13  1.45477E-64     medium             51.82448562  1   6.06902E-13     ethnicity          66.43311272  4   1.28372E-13     gender             9.557732159  1   0.001991089     age                8.182416401  1   0.004229838     enrollment         0.149034707  1   0.699459539     Pell               1.435801885  1   0.230819897     TANF               13.54535541  1   0.000232867     GPA                179.7209374  3   1.01328E-38     onlineExp          12.78010977  3   0.005137086     level                           0        type                            0        motivation         18.00342825  3   0.000439134     medium:level       5.328089832  1   0.020984371     medium:type        5.98948654   1   0.014391391     medium:motivation  6.747183841  3   0.080407666      For the summary output, it looks totally normal, except that none of the values for level or type appear in the list of coefficients.


title: How does the beta prior affect the posterior under a binomial likelihood

text: I have two questions, Question 1: How can I show that the posterior distribution is a beta distribution if the likelihood is binomial and the prior is a beta Question 2: How does choices the prior parameters affect the posterior? Shouldn't they all be the same? Is it possible to answers these questions in R?


title: Application of central limit theorem

text: I'm wondering if I could use an CLT theorem to show that $S_n = \sum_{i=1}^n a_i X_i,$ converges to a normal distribution for large $n$, where $a_i$ are real constants and $X_i$ are I.I.D. poisson random variables. What are the conditions that this approximation is possible? Thanks, Smithy


title: Classification model to predict continuous classes

text: Is it possible to create a classification model which can predict continuous classes, like a number? So far I've working with predictor which can predict one of two classes. I've searched about it and I just found references about continuous features. I have a bunch of features (numbers) to predict a variable, namely 12 features, and observations labelled (with the variable to predict) that can be considered train data.


title: What is the probability of winning this Jackpot Lottery?

text: I am working on some more exam practice questions, and just want to see if I am approaching this question correctly? Here is the complete question: 2,250 numbers are drawn at random without replacement from a barrel containing 150,000 ticket numbers. After all the 2250 numbers are drawn, they are returned to the barrel and a second number (a jackpot number) is drawn at random. If this jackpot number is the same as one of the 2,250 previously drawn numbers, a jackpot prize is awarded. If the jackpot prize is NOT won (that is, the prize winning number is different from one of the 2,250 jackpot numbers), then the number is returned to the barrel and another number is drawn at random. The selection process is repeated until the jackpot prize is won. (1) Find the probability that the jackpot prize is won at any draw. (2) The jackpot prize is initially 10,000, and it increases by 10,000 each time the jackpot prize is NOT won. Find is the probability that the jackpot prize will exceed 250,000 when the jackpot is finally won. (3) Let X be a random variable representing the number of draws until the jackpot prize is finally won. What is P(X = k)? What is the corresponding total prize money won? (4) What is the expected winning of this Jackpot Lottery? Here is my train of thought: (1) This one seemed like a bit of a trick question - since it explicitly states that the selection process is repeated until the jackpot is won - so the probability should be 1. Otherwise, I believe that this can be modeled by a geometric discrete distribution? So the probability of the jackpot prize being won at $x$, should be: $P(X=x)=(1- \frac{2,250}{150,000})^{x-1}(\frac{2,250}{150,000})$ (2) The first "success" in which the prize money exceeds 250,000 would be on the 26th try and thereafter the prize money exceeds this amount, so the probability of the jackpot prize exceeding 250K is: $1-P(X \le 25)$, which works out to be $0.0298$ (3) I have this answer as essentially the same as (1) (Actually what is the difference?) $P(X=k)=(1- \frac{2,250}{150,000})^{k-1}(\frac{2,250}{150,000})$ and the prize money would be $10,000 \cdot k$ (4) For this question, I derive the expected value for a geometric distribution to be: $U_{x} = \frac{1}{p}$ and since $p = \frac{2,250}{150,000}$ the expected draw is $\frac{150,000}{2,250}$ = 66.67. So the expected winning is $10,000\cdot67$ (rounding up) I would greatly appreciate any feedback for the way I have answered the questions - have I approached it correctly? In particular, what would be the difference between Q1 and Q3, and is my answer to Q4 correct? Many thanks in advance!


title: How to replicate a plot of means for a 2 by 3 by 4 design in R?

text: I have a graph that I did in SPSS ![enter image description here](http://i.stack.imgur.com/FGTir.png) I would like to replicate this in R. The data file is here: http://dl.dropbox.com/u/22681355/sendergraphR.csv UPDATE: I've figured out how to to it for the other graph:               graph2<-read.csv(file="SendergraphR.csv")     ddr <- recast(graph2,no.GREEN+no.RED+Senderidentity~variable,         fun.aggregate=mean,id.var=c("no.GREEN","no.RED","Senderidentity"))          qplot(x=no.GREEN,y=Mean_Message,data=ddr,colour=Senderidentity,          group=Senderidentity,geom="line")+facet_wrap(~no.RED,ncol=1)      Now what I would like to do is to separate this second graph into two columns one looking at the cases where the variable urn1 is blue and the other where its red: here's a graph of how it should look: ![enter image description here](http://i.stack.imgur.com/nX4vd.png) I was thinking of doing it the following way:               ddr <- recast(graph2,no.GREEN+no.RED+Senderidentity+urn1~variable,         fun.aggregate=mean,id.var=c("no.GREEN","no.RED","Senderidentity", "urn1"))     qplot(x=no.GREEN,y=Mean_Message,data=ddr,colour=Senderidentity, group=Senderidentity,geom="line")+facet_wrap(~no.RED,ncol=1) +         facet_wrap(~urn1, ncol=2)      But this doesn't seem to work. What am I doing wrong?


title: How to plot a stair steps function with ggplot?

text: I have the graph like this: ![enter image description here](http://i.stack.imgur.com/loNKw.png) R code for generating it is:               DF <- data.frame(date = as.Date(runif(100, 0, 800),origin="2005-01-01"),                       outcome = rbinom(100, 1, 0.1))     DF <- DF[order(DF$DateVariable),] #Sort by date     DF$x <- seq(length=nrow(DF)) #Add case numbers (in order, since sorted)     DF$y <- cumsum(DF$outcome)     library(ggplot2)     ggplot(DF, aes(x,y)) + geom_path() + #Ploting     scale_y_continuous(name= "Number of failures") +     scale_x_continuous(name= "Operations performed")      I want something like this: ![enter image description here](http://i.stack.imgur.com/Fh9X0.png) The difference is in the form of steps in case of failure (I need rectangular). My questions are:   * How to achieve this with ggplot2?   * Is there any better option for visualizing failure rate over time/iterations performed?   * What variant is easier to understand: this or that or, maybe something different?


title: Is it acceptable to use linear regression to predict a 4 point likert scale from a 7 point likert scale?

text: My dependent variable is measured on a 4 point likert scale and independent variable is measured on a 7 point likert scale. Is it appropriate to run regression analysis on such data with varying lengths of likert scale, especially a 4 point likert scale against a 7 point likert scale.


title: Scaling a distribution in JAGS

text: If I have the following model `y~s*dist(a,b)` where `s` is a constant and `dist` is a given distribution (assume that the parameters `a` and `b` cannot be transformed to account for scaling), how do I model this in JAGS? I get the following error: **"[...] is a logical node and cannot be observed"** simply trying the above description. **EDIT:** I am trying to estimate the parameters of a CIR process (it has a scaled non- central chi square distribution, cf here )               model {         for (i in 2:length(timeSerie)){             ncp[i] <- c*et*timeSerie[i-1]             timeSerie[i] ~ dnchisqr(df,ncp[i])/c         }              #parameters         et <- exp(-kappa*timeStep)         c <- 4*kappa/(sigma^2*(1-et))         df <- 4*kappa*mu/sigma^2              #CIR parameter priors          lkappa ~ dunif(-10,-2)         kappa <- exp(lkappa)         mu ~ dunif(3,10)         lsigma ~ dunif(0,5)          sigma <- exp(lsigma)     }


title: Average per line vs. Average per column: is total average the same?

text: _I'm a student in humanities---maths aren't my strong points, but I'm trying hard, so bear with me please:_ I have a survey with seven questions (A to G) and a population divided in four categories (1 to 4). The answer to each question may be True or False. First, I made a calculation of the total of true answers per person and calculated the average of True answer per person per category (i.e. Total of True answers per line for each individual of a given category, divided by the total number of people for that category). Second, I calculated the total of True answers per question and calculated their respective index of True answers (i.e. Total of True answers per column for each individual of a given category divided by the total number of people for that category). I expect the total of the 7 indices per category to provide the average of True answers in that category. But it happens that my two results differ, so either I am really confused or I made a mistake somewhere. Could you help me clear this up? Thanks. I've tried to put together a simplified example:                    A B C D E F G   T      1.  0 1 0 1 0 0 1   3       2.  1 0 0 1 0 1 0   3      3.  1 0 0 1 0 0 0   2      4.  0 0 0 0 1 0 0   1      5.  1 0 1 0 0 0 0   2           T   3 1 1 3 1 1 1  11      Total population: Tp = 5   Total questions: Tq = 7   Total correct answers: Tc = 11   Total of correct answers for question A: TA = 3   Index of correct answers to question A: (iA) = TA / Tp **And here is where the confusion reigns:**   Average of correct answer per person: = Tc / Tp **=?** iA+iB+iC+iD+iE+iF+iG Now that I write it down again, it should be correct, shouldn't it? If so I must have made a mistake on my spreadsheet, which is worth another question in itself.


title: Repeated measures over time with small $n$

text: I was given data to analyze for a study looking at the effects of a treatment on iron levels at four different time points (before treatment, the day treatment ended, 4 weeks after treatment, and 2-4 months after treatment). There is no control group. They are looking to see if there are significant increases in iron levels at each of the 3 post-treatment time points to the before treatment (baseline) level. Eleven patients had baseline levels but only 8 patients had complete data for all 4 time points ($n$=11, 10, 9, and 8 for each time point). Not only were iron levels measured, but two other laboratory measures were taken at each time point to be compared to baseline. I have a few questions as to how to analyze this. I first thought an RM ANOVA would be appropriate to analyze this data, but I was concerned about the small sample size, the loss of data, and the non-normal distribution of the data. I then considered comparing each post-treatment measure to baseline using Wilcoxon signed-rank tests, but then I run into the issue of multiple comparisons. However, I have read some literature that downplays needing to run multiple comparisons. So overall, I'm dealing with small sample sizes, incomplete data, and multiple comparisons (and whether or not its necessary). I hope this all made sense. I'm new to CrossValidated and was directed here by a colleague as a place to learn from experienced statisticians, so I'd appreciate any advice! Thanks! * * * **Edited to add raw data from comment:** There are four total time points and the outcome variable is continuous. For example, the results at each time point look similar to this:                Baseline (n=11): [2, 7, 7, 3, 6, 3, 2, 4, 4, 3, 14]       1st Post (n=10): [167, 200, 45, 132, ., 245, 199, 177, 134, 298, 111]      2nd Post (n=9):  [75, 43, 23, 98, 87, ., 300, ., 118, 202, 156]      3rd Post (n=8):  [23, 34, 98, 112, ., 200, ., 156, 54, 18, .]


title: Fit Inverse Gamma Distribution to dataset in R

text: I am having some problems in modifying Procrastinator's code that assesses the fit of the Inverse Gamma distribution to some randomly generated data. I am refering to the code appearing in the following link: 'Fixing' PearsonFitML to fit to a Pearson V distribution In my case, I have the dataset save in a .csv file which I directly import in R for further processing. MY modification of this code is as follows:               # Required packages     library(MCMCpack)          data=my.csv.data$V1     hist(data)          # log-likelihood     ll = function(par){     if(par[1]>0 & par[2]>0 & par[3]<min(data)) return( -sum(log(dinvgamma(data-    par[3],par[1],par[2]))) )     else return(Inf)     }          # MLE     mle = optim(c(5,2,2),ll)     params = mle$par          # Fit     hist(data,probability=T,ylim=c(0,2.5))     points(seq(2,4.5,0.001),dinvgamma(seq(2,4.5,0.001)-params[3],params[1],params[2]),type="l",col="red")      It seems that I cannot attach the .csv file in here, but I can send it to anyone who is familiar and willing to help. I get the following error:               Error in optim(c(5, 2, 2), ll) :        function cannot be evaluated at initial parameters      and I cannot understand why due to my lack of knowledge in R.


title: Negative value in density deconvolution

text: I have a set of samples from measurement. It can be expressed as random variable $Z = X + Y$. So that $Z$ is observable, but $X,Y$ is unobservable. According to prior knowledge and histogram of $Z$, I can assume that $X$ following exponential distribution with a manual selected parameter. The estimation of probability density function(pdf) of $Y$ is what I wanted. I use fft method to process this density deconvolution problem. More details can be found in this post. To solve this numerical density deconvolution problem we adopt "deconv" function supplied by pracma package in R. The results contain a block of negative value which conflict the definition of probability density. If I ignore all of those invalid value by setting them to zero and convolute it back, there would be a great mount of differences between reconvoluted value and samples. Why? How to process it? Furthermore, if I want to get parameters of $Y$ after its pdf was estimated, perhaps a mixture of normal distribution, what should I do? ================================================= The background is provided as complementary: The samples are collected Round Trip Time(RTT) marked as $D_{rtt}$. It can be expressed as $D_{rtt} = D_{mech} + D_{delay}$, where $D_{mech}$ refers to time delay caused by mechanism, $D_{delay}$ is the real delay which is the main task for this analysis. $D_{mech}$ can be concluded as mixture distribution of exponential distribution and unifrom distribution. Fortunately, there is a gap between the two basic distributions in the pdf. So that I can process it separately as I modeled above. ================================================= I add up my code for sample processing in R below.               delay_est = function(x, samples = 1000, thresh = 0.1){           library(pracma)         lp = subset(x, value < thresh)         d = hist(lp$value, samples)         len = length(d$density)              fit = function(lambda){           di = dexp(d$mids, lambda)           dx = deconv(c(d$density, rep(0, len)), subset(di, di > 0))           rd = dx$q[dx$q<0]           return(sum(abs(rd)))         }              opt = optim(500, fit, method='Brent', lower=0, upper=10000)         print(str(opt))           di = dexp(d$mids, opt$par)               dx = deconv(c(d$density, rep(0, len)), di)         return(list(dx$q, opt))       }


title: Terminology to Use

text: Should you say a factor is associated with a disease or disease risk? For example, which is better to say:   1. Smoking is associated with lung cancer   2. Smoking is associated with lung cancer risk


title: Gaussian process - dimensionality reduction

text: I saw a method for dimensionality reduction for the squared exponential covariance function (not ARD) whereby one uses a $G\times D$ projection matrix $P$ ($G < D$, $D$ = dimension of the inputs) such that the squared distances are calculated in a lower dimensional space. E.g. the kernel function $$K(x, x') = c \times\exp\left(-\frac{1}{2}P(x-x')^TP(x-x')\right)$$ I'm trying to use this in my implementation of Gaussian process regression, but I'm not sure how to get the gradient with respect to $P$ for optimization. I need the derivative of the marginal likelihood w.r.t. $P$ to maximize (minimizing the negative log marginal likelihood in actuality). Has anyone seen this method before or have any suggestions?


title: Procedure for finding boundaries

text: I have a data set of two variables as class and test score. Where class is a classification made by experts and test score ranges 0-10. The data set looks like:               Class   score     1        5.61     1        4.23     2        6.78     3        8.34     2        7.42     3        9.59      How can I find the boundaries in test scores which match classification? I'm looking for something like,say, > "0-5.15= class 1; 5.16-7.81= class 2; 7.82-10= class 3"


title: Joint probabilities for continuous and binomial variables

text: I am having a dreadful time making progress with the following problem, and would greatly appreciate any pointers you had. The problem states that a device's risk of failure within 8 hours of activation, X, varies from day to day according to the probability density function $f(x) = 20(1-x)/9$ for $0.1 < x < 1$, and $f(x) = 20x$ for $0 < x < 0.1$. I need to find both the conditional and the marginal probabilities of exactly 3 failures within 8 hours of activation among the 4 devices activated today. I am to assume that the number of failures follows a binomial distribution. With other problems, we are given f as a function of both x and y, and it is usually straightforward to find the conditional pdf as $h(y|x) = f(x,y) / f_1(x)$ and the marginal distribution as $f_1(x) = \int_{\infty}^{\infty} f(x,y) dy$. I assume I need to do something along these lines, but obtaining $f(x,y)$ from my specific information is giving me a headache. Instead of x in the equation above, I can put ${\rm Binomial}(x,n,p)$ but, since I know $x$ and $n$, I still only have one variable - $p$. I would be very grateful to anyone who could get me started with this, please. Thank you in advance.


title: R equivalent to cluster option when using negative binomial regression

text: I am trying to replicate a colleague's work and am moving the analysis from Stata to R. The models she employs invoke the "cluster" option within the nbreg function to cluster the standard errors. See http://repec.org/usug2007/crse.pdf for a fairly complete description of the what and why of this option My question is how to invoke this same option for negative binomial regression within R? The primary model in our paper is specified in Stata as follows                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)      and I have replaced this with               pday<-glm.nb(cntpd09~logpop08+pcbnkthft07+pccrunion07+urbanpop+pov00+pov002+edu4yr+     black04+hispanic04+respop+as.factor(pdpolicy)+as.factor(maxloan)+rollover+     as.factor(region),data=data[which(data$isser4 != 1),])      which obviously lacks the clustered errors piece. Is it possible to do an exact replication? If so how? If not, what are some reasonable alternatives? Thanks [Edit] As noted in the comments, I was hoping for a solution that didn't take me into the realm of multilevel models. While my training allows me to see that these things should be related, it is more of a leap than I am comfortable taking on my own. As such I kept digging and found this link: http://landroni.wordpress.com/2012/06/02/fama-macbeth-and-cluster-robust-by- firm-and-time-standard-errors-in-r/ that points to some fairly straightforward code to do what I want:               library(lmtest)     pday<-glm.nb(cntpd09~logpop08+pcbnkthft07+pccrunion07+urbanpop+pov00+pov002+edu4yr+      black04+hispanic04+respop+as.factor(pdpolicy)+as.factor(maxloan)+rollover+      as.factor(region),data=data[which(data$isser4 != 1),])     summary(pday)          coeftest(pday, vcov=function(x) vcovHC(x, cluster="state", type="HC1"))      This doesn't replicate the results from the analysis in Stata though, probably because it is designed to work on OLS not negative binomial. So the search goes on. Any pointers on where I am going wrong would be much appreciated


title: Using lmer() with missing levels

text: I believe I need to run a mixed-effects ANCOVA, and I am having difficulty, presumably, because I have missing data in one of my fixed effects or an incomplete design. I have a random grouping variable, `abs.id`, two categorical variables, `leaf.species` and `cond.time`, and one continuous covariate, `day`. `abs.id` includes up to four individuals `leaf.species` has three levels: **map** , **syc** , and **oak**. The problem arises in `cond.time`, which has 2 levels ( **30** and **150** ) for **syc** and **oak** , but only one level ( **30** ) for **map**. I attempt to run the following model with `lmer()` ( _lme4_ ) in R: `lmer(mass ~ leaf.species*cond.time*day + (1| abs.id), data = growth.data)` And receive an error. So, generally, how do I deal with missing data in an `lmer()`?


title: Small dimensional classification (< 20 features), one (or two) dominant predictors

text: I am trying to learn a logistic regression classifier using glm on a dataset of <20 features, but lots of samples. One of the features is a very strong predictor. As a result, the trained model is predicting extreme probabilities of 1.0 and 0.0 on majority of test data. Although the model converges, I get repeated warnings of such extreme predictions. Q1. Are the posterior probabilities/predictions still valid? Q2. How should I reduce the effect of the strong predictor so that the contributions from other variables are taken into account during inference? Right now, they are completely overpowered by the one predictor. Thanks


title: How do I call nonconsecutive rows in an R dataframe?

text: Still learning basic syntax in R. Basic example code is provided below. I know how to call consecutive rows but what if, for example, I wanted the 1st and 3rd row? Or better yet, since that is easily called by w[-2,] in this example, what if the data set was larger and I had the need to investigate the 3rd, 5th and 8th row?               > w<-data.frame(names=c("betty","freddy","sammy"),high.fiving.ability=c(50,50,100))     > w[2:3,]        names high.fiving.ability     2 freddy                  50     3  sammy                 100      BTW I looked at other posts including this one: Filtering a dataframe but what I am looking for is even more basic than that.


title: covariate significant

text: I'm a bit lost at the results part of my thesis. I have conducted an generalized maximum likelihood (GML) univariate test with Y variable: dependent variable one X variable: independent (consists of 2 conditions) I have a lot of covariates to control for. So i performed different ANOVA's with just one covariate at a time. Some covariates turn out to be significant AND the X variable is also significant.   Some covariates turn out to be significant but the X variable is not significant. The problem:   I previously wrote in the results section of my thesis, that when a covariate turns out to be significant and the X was less significant than in the analysis without the covariates, this means I found a main effect but it did not change the interpretation of the results. But what is the case if the some covariates are not significant, but the x variable is significant and even more significant than in the analysis without covariates? What can I conclude from this finding? I tried to google it but i couldn't find it. I hope you can help me. A lot of thanks in advance.


title: Logistic regression shows linear margins for different levels of independent variables

text: Having read about logistic regression, I understand that the probabilities of the success of the DV, `P(Y=1)`, do not necessarily grow the same at every level of the independent variable. This is why you can't interpret the coefficients right away like for OLS, saying an increase in one unit of `x1` increases `P(Y=1)` by `b1` (=the coefficient of `x1`). Using `margins` and `marginsplot` in Stata, however, can help interpret the results. However, I'm getting a **linear relation** in my marginsplot between P(Y=1) and some (maybe every..haven't checked them all) independent variables. What does this tell me? Is this suggesting, that I'm doing something wrong or indicating that I should look into something again. I'm just curious, because usually you would expect margins to be somehow not linear (that's why you do the marginsplot in the first place)..and this perfect linear relationship kind of sets off my alarm bells :)


title: Repeated-measures multiple regression that takes care of sphericity in R

text: I'm looking for a way to run a repeated-measures multiple regression in R, which would take care of sphericity - either by applying some corrections (such as Huynh-Feldt), or by avoiding the problem in some other way. I have 2 factorial repeated measure variables: 3- and 2-level (`roi_ant`, `roi_lat`), and a quantitative between-subject variable (`pred`), and one dependent quantitative variable (`mv`). I want to test for a full model including all possible interactions between the two within-subject and the between-subject variables (i.e., `mv ~ pred * roi_ant * roi_lat`). I am most interested in the slope of `pred` \- whether it is different from 0, and whether it changes depending on `roi_ant` and `roi_lat`. I have tried to do repeated-measures MANCOVA (that would remove the sphericity assumption) using car::Anova package, but if I got this discussion right, `car::Anova()` is not able to handle such a design. `ezANOVA()` performs ANCOVA with applying corrections for sphericity, but reports only the repeated-measure variables (after removing the influence of the covariate), and does not report anything on the covariate `pred`. Mixed linear models (`lme4`) should handle this design well, but since `mcmsamp` is still not implemented, it is difficult to get p-values out of them. Do you have any other suggestions?


title: Need to verify sample is representative of known population of animals

text: I'm a biologist trying to verify that the age, gender, and family group of my samples from a population of whales is representative of the population, in order to make inferences in the discussion. I've been in the field a long time and cannot remember what tests I use to compare goodness of fit in a known population. I am also confused if I should be using two different tests for the continuous variables (age) and the categorized ones (gender and family). My data are below. Could someone give me a clue?                      Samples  Population     Gender  F   18  49             M   11  36             UK  3   3          Family  J   13  26             K   7   20             L   3   42             UK  9   0      ![age distribution here](http://i.stack.imgur.com/6K1VK.png)


title: Diagnosis of normality assumption and fitting linear regression in R

text: I am using the data `veteran` from `R` package `survival`. How can I diagnose the normality assumption about `time`? Should I need to perform a linear regression to measure the dependency of `time` on `age` and `karno`? Are the following commands all to answer the question? or Do I need to add more information?                  #diagnoses normality assumption about 'time'.        shapiro.test(veteran$time)        qqnorm(veteran$time)        qqline(veteran$time)             #perform a linear regression        with(veteran,lm(time~age+karno))        with(veteran,cor.test(time,age))        with(veteran,cor.test(time,karno))      ### EDIT : I want to diagnose the normality assumption. How can I do this for the above data? And for diagnosing the normality assumption, I chose `R` software. So my question involves both diagnosie the normality assumption and R syntax.


title: Distribution of deterministic terms when process is non-stationary (integrated of order I(1), has a unit root)

text: I am looking for literature (or a solution if there is an easy one) on the asymptotic distribution of deterministic terms in linear regression models when the process is unit-root non-stationary I(1). In particular the model I am interested in (set $\alpha=1$ so the process is I(1)): $y_t = \alpha y_{t-1} + \delta step_t + \epsilon_t$ where step_t is a step function for a couple of periods. Most of the literature I have found focuses on the distribution of $\alpha$ and testing whether it is equal to $1$ or not. What I am trying to figure out is the distribution of $\delta$. It will be asymptotically normal if $\alpha < 1$, but what if \alpha=1? So overall the question is, is there a literature out there that addresses the asymptotic distribution of deterministic terms when the process is unit-root non-stationary? Thanks!


title: How To Write Seasonal ARIMA model mathematically

text: I'm trying to write a seasonal ARIMA model ARIMA(1,0,3)(1,2,0) period 5 mathematically but I don't seem to be able to follow what this resource is saying otexts arima The example they use is ARIMA(1,1,1)(1,1,1)4 so it makes it very hard for me to follow what's actually going on because all the numbers are the same. Could someone help me write out my model mathematically and possibly explain how it's done?


title: Question about Time 1/Time 2 analysis

text: I was hoping to get some input on the best way to test some data that I have collected (using SPSS). In short, respondents make an evaluation at time 1, then complete a distractor task, and then complete an evaluation at time 2 (DV is a 1-7 scale). It is a 2x2 between subjects design where 2 variables are manipulated at time 2. There are likely several ways to approach this, but do you all have any suggestions for what may work best to analyze the differences between means at time 1 and time 2? And I suppose it is technically a repeated measures design, as respondents make 2 evaluations each. Any thoughts would be very much appreciated! Thanks...


title: Which statistical method should be used for finding statistical significance of any observed correlation?

text: My knowledge of statistical methods is so poor, so I don't know exactly what is the best statistical method to opt in my case. **Ultimate goal: These parameters(A, B and C) could affect parameter D. I want to know whether they do affect it? What is the correlation between each of the parameters and parameter D and how significant this correlation is.** Imagine, I measured D (7 times during a year). Each time, I also measured A, B and C quantities of each sample (I have three different samples). So now, I have C (dependent variable) vs A, B and C parameters. What I want to know is the statistical significance of any observed correlation between the dependent variable and the independent variables. My null hypothesis to be tested is 'there is no correlation between the dependent and independent variables.' And I want to investigate these significances by means of defining three confidence level of 90%, 95% and 99%. example of expected output: the correlation between A and D is Significant for 99% and 95% confidence level but not for 90% confidence level. Could anyone help me how to perform this in Excel?


title: ordinal regression or Spearman correlation?

text: I am conducting research investigating two different groups of women: one group of Short-Stature Women (SSW) and another with Non-Short Women (NSW). We have the hypothesis that SSW has an inaccurate auto-perception of their current body size (CBS). We assessed CBS with a figure rating scale, consisting of 9 different silhouettes, ranging from 1 to 9 (which would assess something like from "malnourished" to "very obese"). We know that CBS is highly determined by the person's Body Mass Index (BMI). Hence, we would like to know if the BMI is a significantly better "predictor" of CBS in NSW compared to SSW. My first approach was to run two Spearman rho rank correlation test (one for each group (SSW and NSW), between BMI and CBS) and then to compare both coefficients using the Fisher's Z-test. Nevertheless, I am pretty sure that I can treat CBS as an ordinal variable, right? So, it seems to me like an ordinal regression (using CBS as DV and BMI as IV) would be more suitable to my data. If this is true, my problem is that I do not know how to compare coefficients for ordinal regression. Would it be comparing Nagelkerke-R²? How can I do this? To summarize I have 2 main questions: a) Is Spearman rho rank correlation with comparison of the coefficients by Fisher's Z-test adequate in this case? and b) If ordinal regression is more adequate, how can I compare the coefficients of the regression between my two groups?


title: Measure of closeness between samples (given a larger sample in the background)

text: I have data in a 15d-space and want to evaluate closeness between samples of instances. As i dont have any assumptions on the distribution of the data I prefer Bhattacharyya coefficients over Mahanolobis distance. But as my samples are pretty small, I am unsure whether Bhattacharyya coeficients will return reasonable distance estimates. I'd like to compare samples of size ~100vs.~4 and sometimes even ~4vs.~4. As the pool of data I am drawing from is pretty large, I was wondering if there is a method that is using this additional knowledge I have about the distribution, too. Even if I have to recover the distribution from the additional data by some kernel density estimation I would be happy. If there is some parameter-free method out there it would be even better. Kind regards Damian


title: Nonparametric regression with restrictions

text: I need to estimate the function $f$ in the following object $$E[Y|X=x] = a(x)f(x) + b(x)f'(x)$$ with $a(x),b(x)$ known. What would be the optimal strategy to do this? I tried to do the cubic spline approximation with two knots: $f(x) = \theta_0+\theta_1x+\theta_2x^2+\theta_3x^3+\delta_1(x-\xi_1)_+^3+\delta_2(x-\xi_2)_+^3$ and so $$E[Y|X=x] = \theta_0a(x)+\theta_1(a(x)*x+b(x))+\theta_2(a(x)*x^2+b(x)*2x) + \theta_3(a(x)*x^3+b(x)*3x^2)+\delta_1[a(x)*(x-\xi_1)_+^3+3*b(x)(x - \xi_1)_+^2] + \delta_2[a(x)*(x-\xi_2)_+^3+3*b(x)(x - \xi_2)_+^2]$$ where all the parameters are estimated with OLS. But it does not work numerically. Increasing the number of knots and choosing a different type of series approximation, e.g. B-splines does not help. What am I doing wrong?


title: Grouping of procedure codes and diagnosis codes

text: I am trying to group procedure codes so that I do not need to hard code thousands of procedure codes and diagnosis codes. I am trying to build a regression to see the effect of percentage of dollars for a provider on total billed amount for a provider on several parameters so that I can find the outliers. The parameters are procedure codes, diagnosis codes, place of service etc. Any suggestions on grouping procedure codes and diagnosis to see the effect on total billed amount for a provider would be highly appreciated.


title: Accuracy vs. simplicity as criteria to select explanatory variables

text: I want to create a binary model which predicts whether someone has improved his state. I am testing possible variables as explanatory variables in order to make some recommendations. Now my binary model is quite weak, AUC is around 0.6 and there are hardly any cases predicted above 0.5 in the test set. Now I could simply improve my model by including the previous state which is highly correlated with improvement, but for me this is too obvious. It might 'kick out' other explanatory (more actionable) variables. So now I am faced with this question of whether to go for accuracy or simpler/actionable model. What do you recommend?


title: Cluster size in generalized estimation equation (GEE)

text: I'm using the `gee()` function from the `gee` package in R. The problem I'm having is that the 'Maximum cluster size' that I get from the output of the GEE function seems to disagree with what I believe it should be given my data. Here's a small example, where I have six observations from each of ten patients, which are members of one of three groups:               ID <- rep(1:10, 6)     myData <- data.frame(ID)     myData$ID <- factor(myData$ID)     myData$Group <- c(rep('A', 20), rep('B', 20), rep('C', 20))         myData$value <- rnorm(60, mean = 10, sd = 1)      Calling `summary(myData$ID)` shows that I do, in fact, have six observations of each ID, and calling `class(myData$ID)` shows that `ID` is, in fact, a factor. Therefore, I would expect that the maximum cluster size would be six. However, when I call the following...               gee(value ~ Group, id = ID, data = myData)      ...the `Maximum cluster size` that prints out is `1`. Am I misunderstanding what maximum cluster size means? Incorrectly formatting my data? Thanks in advance.


title: Error messages when doing regression in SPSS

text: I need to do a number of regressions in SPSS but keep getting error messages. For example - "There are X cells (i.e., dependent variable levels by combinations of predictor variable values) with zero frequencies." and "The log-likelihood value is practically zero. There may be a complete separation in the data. The maximum likelihood estimates do not exist." Example of one of the regressions I am trying to carry out - I'm examining carcass damage in pigs, specifically, the effect of injuries that occurred on farm on the injuries visible on the carcass. I want to predict skin lesion severity at the abattoir (absent, mild, moderate, severe) from skin lesion severity on farm (absent, mild, moderate, severe) and the number of weeks the animals had moderate or severe skin lesions (0 weeks, 1 week, 2 weeks, 3 weeks). I want to do an ordinal regression on this data. Could somebody explain to me in simple terms how to avoid the "zero frequencies" error message and other error messages? Update - I have researched this online and now understand why I am getting the "zero frequencies" message. If there are some zero frequencies does this mean that I can't do a regression on the data? Is there an alternative?


title: Completing a matrix whose entries depend on time

text: Imagine you have a matrix, where each row represents an individual and each column represents a specific time interval. The entries of the matrix represent events occurring for a given person at a specific interval of time. Now say we wish to use this representation to build a classifier. So along with the matrix described above, we have a vector of labels for an event of interest, such as death. Now, how would one deal with making predictions for this representation on individuals who have been observed for varying lengths of time? That is, how do we fill the entries of the matrix, given that we cannot have observed all of the possible time intervals for each instance?


title: Association rules for continuous features

text: I need to detect rules in a data set that contains real-valued features. A simplified example of my samples, defined by features a, b, and c, and having class 0 or 1, might look like this: > 2.34, 1, 5.46 => 0 > > 1.23, 0, 7.81 => 1 For this data set, I am interested in detecting rules like:               a > 1 and a<3 and b==1 => class=0      I have checked out association rule algorithms in Weka and then searched online, and they do what I want, but they all seem to need binary features. My question: is there an algorithm that can do what I want, or do I need to hack my own on the basis of existing association rule algorithms?


title: Initial values for logistic regression using maximum likelihood

text: I'm trying to calculate logistic regression coefficients by defining the log- likelihood function and using maximum likelihood. In some cases when the initial (start) values I gave to the maximum likelihood were not correct I got wrong results for the logistic regression (different from the ones I get when using `glm` for example). Given the input data and y values, what should be the optimum initial values for logistic regression (or, in other words, what are the values that are being used in `glm`)?


title: Prior Gamma distribution: Select appropriate alpha given beta and median

text: I am trying to programatically select a prior distribution from the Gamma family of distributions. The primary criteria that I need to satisfy is that the median of the distribution should be a given value $X$ (i.e. such that it's equally probable to have a parameter value above $X$ or below $X$) Additionally, because I'm using the distribution as a conjugate prior for a Poisson distribution's lambda, and I know how much relative weight I want to give to the prior relative to the subsequent observations, I _think_ I know what an appropriate $\beta$ parameter is (though I could be wrong on this point, feel free to make comments if I'm thinking about that incorrectly). So essentially my problem becomes how to select an appropriate $\alpha$ such that the median is $X$. The way I've considered this seems to work well to get a good prior distribution when constraining the mean (or even mode), but technically, I should be constraining the median. Is this possible formulaicly, or would I need to resort to some kind of computational search to find the right value?


title: Forecasting Interventions (pulse) with ARIMA Model

text: **Q1:** Is there any Arima (p,d,q) model that can forecast interventions (pulse) itself? I know that I can use `xreg` or even `xtransf` arguments as the covariates to include the intervention over the observed time series. The problem is that I don’t know the new values of the covariates to include in the predict function to do the forecasts. **Q2:** Can we simulate (estimate) these new covariates? If yes, how? I think that depends on the probability of having pulse. But how we can estimate these probabilities? To be more precise, I want to develop a time series model that can capture the interventions and also be able to forecast them! Is that even possible? Maybe I am asking too much :) I appreciate any comment, suggestion or reference.


title: Composite Scores and Standardized Composite Scores t test

text: I have a set of survey data related to 20 survey questions. Each of these questions represent a variable (`Q1`, `Q2`,...`Q20`). I created a new variable `QCom` which measures the response of the survey, and is given by a composite score obtained as the sum of the scores of responses from `Q1` to `Q20`. I then perform a t-test for `QCom` to check for evidence of a difference in mean due to another variable (`Sex` -> Male or Female). I obtained the test statistics from SPSS's independent sample t test. Subsequently, I created another variable called `StCom` (equal to standardized score of `QCom`). Again, I repeated the t-test as the above using `StCom` this time. The t test statistics that I obtained from SPSS was exactly the same as the first test using `QCom`. In this case, I am not sure if this is normal or the z-scores transformation is incorrect. Can someone help to enlighten me why the composite scores and composite z-scores t test results are exactly the same?


title: Network graph: hypothesis test that degree distribution is becoming more equitable

text: I'm playing around with network graphs and am wondering how to show that there is:   1. a statistically significant difference in in-degree distribution between nodes for 2 time periods   2. the in-degree distribution is more equally dispersed among nodes at time 2 Now let's say I have a classroom of students and I asked them to name the 3 kids who they thought was getting the most negative attention from other students. Let's say three kids (J, K, and L) were getting the most negative attention (we'll say they're being bullied). Now we conduct a session (treatment) and now conduct a follow up question to the classroom again asking the three students that were receiving the most negative attention. Here I plotted that out with R as a network graph. You can see from the weighting I obviously made this the 'ideal' situation that I have described and after the treatment the vertices are more equally distributed. How can I show a statistically significant: (a) change in the distribution of in-degrees (incoming edges) and (b) the distribution of in-degrees (incoming edges) is more equal? ![enter image description here](http://i.stack.imgur.com/rEMTL.png)               library(igraph); library(reshape2)          weights_1 <- c(.01, .01, .01, .01, .01, .05, .05, .05, .05, .25, .25, .25)     weights_2 <- rep(1/11, 12)          selects <- function(p, n = 12) {         melt(data.frame(from = LETTERS[1:n],                 do.call(rbind, lapply(1:n, function(i) {                     sample(LETTERS[1:n][-i], 3, FALSE, prob = p[-i])                 }))),              id="from", value.name="to",              variable.name = "choice.n"         )[, c(1, 3, 2)]     }          set.seed(100)     time_1 <- selects(weights_1)     t1 <- graph.data.frame(time_1)          set.seed(100)     time_2 <- selects(weights_2)     t2 <- graph.data.frame(time_2)          E(t1)$arrow.size <- .6         E(t2)$arrow.size <- .6          par(mfrow=c(1,2))     plot(t1, layout=layout.circle); mtext("Time 1")     plot(t2, layout=layout.circle); mtext("Time 2")      **EDIT** It occurs to me that it may be easier to help if the Boolean matrices and adjacency matrices are already created: **Here it is as a viewable format:**               $Boolean_time1       A B C D E F G H I J K L     A 0 0 0 0 0 0 0 0 0 0 0 0     B 0 0 0 0 0 0 0 0 0 0 0 0     C 0 0 0 0 0 0 0 0 0 0 0 0     D 0 0 0 0 0 0 0 0 0 0 0 0     E 0 0 0 0 0 0 0 0 0 0 1 1     F 0 0 1 0 0 0 0 0 0 0 0 1     G 1 0 0 0 0 0 0 0 0 1 0 1     H 0 0 0 0 1 0 0 0 1 0 0 0     I 0 0 0 1 0 0 0 1 0 0 0 0     J 1 1 0 1 0 1 1 1 0 0 1 0     K 0 1 1 1 1 1 1 1 1 1 0 0     L 1 1 1 0 1 1 1 0 1 1 1 0          $Boolean_time2       A B C D E F G H I J K L     A 0 0 0 0 0 0 0 0 0 0 1 1     B 0 0 0 0 0 0 0 0 0 0 0 0     C 0 1 0 1 0 0 0 0 1 0 0 0     D 0 0 0 0 0 1 0 0 0 1 0 0     E 1 0 0 0 0 0 1 0 0 0 1 0     F 1 0 1 0 1 0 0 0 1 0 0 0     G 0 0 0 0 1 1 0 1 0 1 1 0     H 0 1 1 0 0 0 1 0 0 0 0 1     I 1 1 0 0 0 0 0 0 0 0 0 0     J 0 0 0 1 0 1 1 1 0 0 0 0     K 0 0 1 0 1 0 0 1 1 0 0 1     L 0 0 0 1 0 0 0 0 0 1 0 0          $adjmat_time1       A B C D E F G H I J K L     A 0 0 0 0 0 0 0 0 0 0 0 0     B 0 0 0 0 0 0 0 0 0 0 0 0     C 0 0 0 0 0 0 0 0 0 0 0 0     D 0 0 0 0 0 0 0 0 0 0 0 0     E 0 0 0 0 2 1 1 0 0 1 0 1     F 0 0 0 0 1 2 1 0 0 0 1 1     G 0 0 0 0 1 1 3 0 0 1 1 2     H 0 0 0 0 0 0 0 2 0 0 2 2     I 0 0 0 0 0 0 0 0 2 2 2 0     J 0 0 0 0 1 0 1 0 2 7 5 5     K 0 0 0 0 0 1 1 2 2 5 9 7     L 0 0 0 0 1 1 2 2 0 5 7 9          $adjmat_time2       A B C D E F G H I J K L     A 2 0 0 0 1 0 1 1 0 0 1 0     B 0 0 0 0 0 0 0 0 0 0 0 0     C 0 0 3 0 0 1 0 1 1 1 1 1     D 0 0 0 2 0 0 2 0 0 1 0 1     E 1 0 0 0 3 1 1 1 1 1 0 0     F 0 0 1 0 1 4 1 1 1 0 3 0     G 1 0 0 2 1 1 5 0 0 2 2 1     H 1 0 1 0 1 1 0 4 1 1 2 0     I 0 0 1 0 1 1 0 1 2 0 0 0     J 0 0 1 1 1 0 2 1 0 4 1 1     K 1 0 1 0 0 3 2 2 0 1 5 0     L 0 0 1 1 0 0 1 0 0 1 0 2      **As an easy to grab R format:**                dats <- structure(list(Boolean_time1 = structure(list(A = c(0L, 0L, 0L,                           0L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 1L), B = c(0L, 0L, 0L, 0L, 0L,                                0L, 0L, 0L, 0L, 1L, 1L, 1L), C = c(0L, 0L, 0L, 0L, 0L, 1L, 0L,                                0L, 0L, 0L, 1L, 1L), D = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L,                                1L, 1L, 0L), E = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L,                                1L), F = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 1L), G = c(0L,                         0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 1L), H = c(0L, 0L, 0L,                                0L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 0L), I = c(0L, 0L, 0L, 0L, 0L,                                0L, 0L, 1L, 0L, 0L, 1L, 1L), J = c(0L, 0L, 0L, 0L, 0L, 0L, 1L,                                0L, 0L, 0L, 1L, 1L), K = c(0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L,                                1L, 0L, 1L), L = c(0L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 0L,                                0L)), .Names = c("A", "B", "C", "D", "E", "F", "G", "H", "I",                                 "J", "K", "L"), row.names = c("A", "B", "C", "D", "E", "F", "G",                              "H", "I", "J", "K", "L"), class = "data.frame"), Boolean_time2 = structure(list(                  A = c(0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 1L, 0L, 0L, 0L), B = c(0L,                              0L, 1L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L), C = c(0L, 0L,                                    0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 0L), D = c(0L, 0L, 1L,                                    0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L), E = c(0L, 0L, 0L, 0L,                                    0L, 1L, 1L, 0L, 0L, 0L, 1L, 0L), F = c(0L, 0L, 0L, 1L, 0L,                                    0L, 1L, 0L, 0L, 1L, 0L, 0L), G = c(0L, 0L, 0L, 0L, 1L, 0L,                                    0L, 1L, 0L, 1L, 0L, 0L), H = c(0L, 0L, 0L, 0L, 0L, 0L, 1L,                                    0L, 0L, 1L, 1L, 0L), I = c(0L, 0L, 1L, 0L, 0L, 1L, 0L, 0L,                                    0L, 0L, 1L, 0L), J = c(0L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 0L,                                    0L, 0L, 1L), K = c(1L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L,                                    0L, 0L), L = c(1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L,                                    0L)), .Names = c("A", "B", "C", "D", "E", "F", "G", "H",                                  "I", "J", "K", "L"), row.names = c("A", "B", "C", "D", "E", "F",                              "G", "H", "I", "J", "K", "L"), class = "data.frame"), adjmat_time1 = structure(c(0,           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,                                0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,                                0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 1, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0,                                1, 2, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 3, 0, 0, 1, 1, 2, 0,                                0, 0, 0, 0, 0, 0, 2, 0, 0, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2,                                2, 0, 0, 0, 0, 0, 1, 0, 1, 0, 2, 7, 5, 5, 0, 0, 0, 0, 0, 1, 1,                                2, 2, 5, 9, 7, 0, 0, 0, 0, 1, 1, 2, 2, 0, 5, 7, 9), .Dim = c(12L,                             12L), .Dimnames = list(c("A", "B", "C", "D", "E", "F", "G", "H",                              "I", "J", "K", "L"), c("A", "B", "C", "D", "E", "F", "G", "H",                                "I", "J", "K", "L"))), adjmat_time2 = structure(c(2, 0, 0, 0,                                 1, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,                                0, 3, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 2, 0, 0, 2, 0, 0, 1,                                0, 1, 1, 0, 0, 0, 3, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 0, 1, 4, 1,                                1, 1, 0, 3, 0, 1, 0, 0, 2, 1, 1, 5, 0, 0, 2, 2, 1, 1, 0, 1, 0,                                1, 1, 0, 4, 1, 1, 2, 0, 0, 0, 1, 0, 1, 1, 0, 1, 2, 0, 0, 0, 0,                                0, 1, 1, 1, 0, 2, 1, 0, 4, 1, 1, 1, 0, 1, 0, 0, 3, 2, 2, 0, 1,                                5, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 1, 0, 2), .Dim = c(12L, 12L), .Dimnames = list(                  c("A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K",                                      "L"), c("A", "B", "C", "D", "E", "F", "G", "H", "I", "J",                                     "K", "L")))), .Names = c("Boolean_time1", "Boolean_time2",                                "adjmat_time1", "adjmat_time2"))


title: updating posterior parameters when involving conditioning

text: I have a setup where the joint posterior is written as: $$ P(w, \lambda, \phi \vert y) = P(\phi) \times P(w \vert \lambda) \times P(\lambda) \times \prod_{i=1}^{N}P(y_i \vert w_i, \phi, \lambda) $$ Now $\phi$ and $\lambda$ are modelled as Gamma distributions and the likelihood and prior on $w$ is normally distributed. So, taking advantage of the conjugate prior property, the posterior on $\lambda$ and $\phi$ should be Gamma distributed. So, my strategy is that for a given approximation for $w$, I will update the posterior parameters for $\lambda$ and $\phi$. So, for example, if I want to update $\lambda$, I would do something like: $$ P(\lambda) = \int \int P(w, \lambda, \phi \vert y) \:dw \: d\phi $$ So, I marginalise out $w$ and $\phi$, evaluate the integral and match the terms to get the desired Gamma distribution parameters. Now, I want to do it iteratively, so for a given instance 't' I have $w=w_t$. So, I would want to update the parameters of $\lambda$ and $\phi$ by conditioning on $w=w_t$. I have two very basic questions about this: 1: In the integral, can I basically replace the terms $w$ with $w_t$ and treat them like constants wrt to the variable of interest? I am slightly weak on the calculus side. This should be the case but I wanted to check. So, $P(\lambda|w=w_t) \propto \int \int P(w_t, \lambda, \phi \vert y) \: d\phi$ 2: Also, $P(\lambda \vert w=w_t) = \frac{P(\lambda, w)}{P(w=w_t)}$ which is equal to $\frac{\int P(w, \lambda, \phi) \: d\phi}{P(w=w_t)}$. However, since the posterior $P(\lambda)$ is Gamma distributed, I am guessing I can always compute the normalisation factor at the end and do not have to worry about the denominator. Is that correct? So, all I want to do at the moment, is update the parameters and compute the expectation of the variable $\lambda$ and I guess the scaling should not matter in this case.


title: Autocorrelation, Durbin-Watson and non time-series data

text: I have a simple linear regression with age as independent variable and a cognitive scale as dependent variable. Each subject is present only once. As it is not time-series data and there is no spatial effect, is it correct not to check for autocorrelation? Does a Durbin-Watson result of .23 mean something?


title: Program to compute partial derivatives

text: I need to find partial derivatives for: $L=\frac{1}{2}\sum_{i=1}^{n} w_{i}^{2}\sigma_{i}^{2} -\lambda \left( \sum_{i=1}^{n} w_{i} \bar{r_{i}} - \bar{r} \right) -\mu \left( \sum_{i=1}^{n} w_{i}-1 \right)$ with 5 variables, I get 5 partial derivatives $\frac{\partial L}{\partial w_{1}}$, $\frac{\partial L}{\partial w_{2}}$, $\frac{\partial L}{\partial w_{3}}$, $\frac{\partial L}{\partial \lambda}$ and $\frac{\partial L}{\partial \mu}$ with which I need to solve for $w_{1}$, $w_{2}$ and $w_{3}$. For all $i$, $\sigma_{i}$ is known. I don't want to manually differentiate; it could be more cumbersome with more restrictions. Later, I need to find the point where the minimum variance is located; i.e. to find out the optimal values for $\lambda$, $\mu$ and $w_{i}$ for all $i$ where $i\in\\{1,2,3\\}$. I need to vary the 5 variables to find the optimal solutions (a bit like in Excel but no such thing currently available) and then show the frontier graphically. So my question is which program would you use to solve such thing?


title: Determine which method predicts best given results from a simulation exercise

text: This question is related to one I have already asked but the answer I got that suggests I should adopt a new tack to address my research question. I repeat the substantive part of the original question to show the features of my particular computer simulation/experiment: > I conducted a computer-based assessment of different methods of fitting a > particular type of model used in the palaeo sciences. I had a large-ish > training set and so I randomly (stratified random sampling) set aside a test > set. I fitted $m$ different methods to the training set samples and using > the $m$ resulting models I predicted the response for the test set samples > and computed a RMSEP over the samples in the test set. This is a single > **run**. > > I then repeated this process a large number of times, each time I chose a > different training set by randomly sampling a new test set. > > Having done this I want to investigate if any of the $m$ methods has better > or worse RMSEP performance. I also would like to do multiple comparisons of > the pair-wise methods. Say I did 50 Runs so I have $50 \times m$ observations of the RMSEP. I wish to determine if all $m$ methods yield the same RMSEP (test a null of equal RMSEP for each method). Assuming a difference in RMSEP between models I would also like to know which models differ in the predictive performance. How might I go about addressing this question statistically? * * * Note that the reason I did 50 Runs was to avoid the issue of getting a particular result just because I chose that particular training set / test set combination. Computationally I can probably afford to do and order or magnitude more Runs quite easily if required.


title: Question about combinations involving a caveat/exclusion

text: Working through an example in the Ross Probability textbook. Can someone explain the reasoning behind the answer? From a group of 7 men, how many committees consisting of 3 men can be formed? what if 2 men are feuding and refuse to serve on the committee together? **Answer, part 1:** First part is clear for me. with n=7 and r=3. Result = $\binom{7}{3}=35$. **Answer, part 2:** Textbook gives the answer as 35 - 5 = 30, with 5 a result of $\binom{2}{2}\binom{5}{1}$. Can someone please explain the reasoning behind the Part 2 answer?


title: Coding single word responses into consistent

text: this might be a stupid question, but I'll try anyway! I have a data set from a survey asking about what brands they can remember within the category toys. The survey participants get to write a single brand in 10 different text boxes. The purpose is to reveal the brand with the greatest top of mind effect, without displaying any suggestions. **The problem is:** text answers are not consistent. For example "Toy's are us" vs "Toys r us'" **My question is:** is there any plug in for Stata correcting similar answers based on a scope of answers? (answers out of this scope can be set to missing) All tips and answers are appreciated!


title: Adaptive regression splines in earth package R

text: I am using earth package for the following data.               x <- c(127, 128, 255, 256, 511, 512, 600, 700, 800, 900, 1000, 1023, 1100,            1200, 1300, 1400, 1500, 1600, 2047, 2048, 2100, 2200, 2300, 2400, 2500,            2600, 2700, 2800, 3000, 3100, 3200, 3300, 3500, 4063, 4064, 4100, 4200,            5200, 5400)          y <- c(0.59, 0.61, 0.59, 1.55, 1.33, 3.50, 1.00, 1.22, 2.50, 3.00, 3.79,            3.98, 4.33, 4.45, 4.59, 4.72, 4.82, 4.90, 4.96, 7.92, 5.01, 5.01,            4.94, 5.05, 5.04, 5.03, 5.06, 5.10, 5.04, 5.06, 7.77, 5.07, 5.08,            5.08, 5.12, 5.12, 5.08, 5.17, 5.18)       After building the model,               model<-earth(y~x)      I get following regression model.               summary(model)     Call: earth(x=x, y=y)     coefficients     (Intercept)  5.225822553     h(1400-x)   -0.003820087      Is there any possibility that I can increase somehow the number of knots or regression splines?


title: Synthetic time series generation

text: I have a linear model (with seasonal dummy variables) that produces monthly forecasts. I'm using R together with the 'forecast' package:               require(forecast)     model = tslm(waterflow ~ rainfall + season, data = model.df, lambda = lambda)     forec = forecast(model, newdata = rainfall.df, lambda = lambda)      I did a cross-validation and it looks great. Now, what i need is to generate _weekly data points_ from these month forecasts - in other words, i need to generate a synthetic time-series that have monthly means equal to the forecasts above. So my function would look like:               generate.data = function(monthly.means, start.date, end.date)     {        #code here     }      I'm not sure how to do this (interpolation?), so any help is welcome. Thanks!


title: Whither bootstrapping - can someone provide a simple explanation to get me started?

text: Despite several attempts at reading about bootstrapping, I seem to always hit a brick wall. I wonder if anyone can give a reasonably non-technical definition of bootstrapping? I know it is not possible in this forum to provide enough detail to enable me to fully understand it, but a gentle push in the right direction with the main goal and mechanism of bootstrapping would be _much_ appreciated! Thanks.


title: Understanding the p-value

text: I'm studying the second course in statistics, and now I have a problem to understand $P$-values. Namely, one exercise is the following. When testing hypothesis $H_0:\mu=\mu_0$ one gets the value of test statistic $z=1.7$. Determine the $P$-value for the alternative hypothesis $H_A:\mu>\mu_0$. Is it that I have to calculate the $P$-value from $N(0,1)$-distribution or is there some other distribution in the background?


title: How to calculate a probability of a sequence of observations in Hidden Markov Model?

text: ![an example of HMM](http://i.stack.imgur.com/fct3Z.jpg) Given a model (from Ankur Jain's slides) and question with the answer is presented below: ![enter image description here](http://i.stack.imgur.com/lX01f.jpg) Initial probabilities: say P(‘Low’)=0.4 , P(‘High’)=0.6 . However i'm confused why in the solution we have an extra 0.4, i guess it should be as follows, can you clarify me?               P(‘Dry’|’Low’)P(‘Rain’|’Low’) P(‘Low’)P(‘Low’|’Low)     = 0.4*0.6*0.4*0.3


title: How to calculate this formula for variance?

text: I have a function which I would like to use Taylor expansion and calculate its variance by the following formula: The formula for variance then becomes \begin{align} \operatorname{Var}(f(X))=[f'(EX)]^2\operatorname{Var}(X)+\frac{[f''(EX)]^2}{4}\operatorname{Var}^2(X)+\tilde{T}_3 \end{align} got the formula from Variance of a function of one random variable I have tried a few times to calculate a simple example to get to learn how to use it, but can not say I have succeeded. So if anyone please could show me how to use the above formula for calculating the variance for this (simple) function: $f(X)=100\times \exp(x)+100\times \exp(2x)$, where the expected value of $x$ is $0.05$ and standard deviation is $0.1$. normal distribution of x is assumed.


title: How to use factors (from CFA) as independent variable in Regression Analysis

text: I calculated 4 factors as latent constructs in a Confirmatory Factor Analysis (I use AMOS). Now I am wondering if it is possible to extract some kind of a factor score like I know it from Exploratory Factor Analysis with SPSS to use those as independent variables in a multiple regression analysis. Can I use the "Data Imputation" Function of AMOS? What are the risks? THank you a lot in advance


title: Modeling Counts With Small Observations

text: I am new to Cross Validated SE so I am going to try and formulate my question to the best of my ability. I have a large data set that contains $5$ different fields. The fields are               Gender-(Male, Female),     Age-(0,18],(18,25],(25,45], (45,65], (65,100],     Region - (Northeast, Midwest, South, West)     MSA - (MSA,Non-MSA)     Doctors Office Visits - Number of times an individual visited the doctors office.       I am trying to make a model of the counts for each possible subset. So for example I would have the subset -               Age - (0,18]     Gender - Male     Region - West     MSA - Non-MSA     Doctors Office Visits - Number of times an individual in this subset visited the doctors office.      Some of these subsets have a good number of observations in which I can model the counts. However, some of the subsets have a small number of observations `(1-3)` which really skews the counts, especially if they are high (or small). I also want to add that there are no missing data in the Doctor Office Visits field. So my questions are:   1. I am using a poisson distribution to model the counts. Is this an effective model? Are there other ways I should model the counts? My goal is to get a predicted number of counts for each subset.   2. For the subsets where there are small observations, what can I do to model those counts so they are not skewed? Are there methods for adding in random data and modeling on that? Or can I use observations from other subsets or possibly the whole data set? I'm not sure how to proceed with this. Any help and comments would be greatly appreciated! Thank you!


title: How do I identify those specific observations in which a given variable has the widest variation in response?

text: My data has at least 5 variables with a significant response, with one variable having the highest regression co-efficient. I want to find those observations where that variable did _not_ behave as expected, to better understand the other contributors to the outcome.


title: Doubt on interpretation of Neyma-Pearson lemma

text: In my book: $\mathbf{X}=(X_1,\ldots,X_n)$ $f(\mathbf{x})$ is the joint density, where $f$ is either $f_0 \text{ or } f_1$. Suppose we want to test $H_0: f=f_0$ or $H_1: f=f_1$. The test, whose test function is $$\phi(\mathbf{X})=1\text{ if }\frac{f_1}{f_0}\geq k;$$ $$\phi(\mathbf{X})=0 \text{ otherwise,}$$ (for some $0<k<\infty$) is a most powerful test of $H_0$ versus $H_1$ at level $E_0(\phi(\mathbf{X}))$. My question is how is $k$ defined? Can I interpret the lemma as if $\forall k \in (0,\infty)$ there will be a test function $\phi(\mathbf{X})$ such that it will determine the size for which the test with test function $\phi(\mathbf{X})$ is most powerful? I'm just trying to understand which kind of relationship $k$ and the size of the test have between each other.


title: how to calculate survival difference between groups at various time points

text: I have two groups of patients which have been followed for up to 5 years. Using SPSS I would like to find the cumulative survival at 1,2 or 3 years and see if there is a statistically significant difference between the groups at each time point. This is the equivalent of right censoring at 1,2 or 3 years and running a KM analysis for each. I tried using life-tables but the Wilcoxon (Gehan) statistic is always the same irrespective of whether I limit the time interval to 1,2 or 3 years.               SURVIVAL TABLE=calc_time_OS BY highLC (1 2)       /INTERVAL=THRU 36 BY 12       /STATUS=SURVIVAL_STATUS(1)       /PRINT=TABLE       /PLOTS (SURVIVAL)=calc_time_OS BY highLC       /COMPARE=calc_time_OS BY highLC.      I am sure there is something pretty straight forward that I am missing. Thanks


title: 3D Gaussian basis function

text: I need tο use 3D data for training and for that I need to find the Gaussian basis function. I know how to find $f(x,y)$ but how can I find $f(x,y,z)$?


title: Alternative to forecast() and ets() in Python?

text: I'm looking for a Python alternative to R's ETS() from forecast(). It's my understanding that ETS() is one of the best performing forecasting program and I would like to use it. However I am extremely uncomfortable using R and I already have all my data cleaned and set up in Python. Does anyone know any other packages similar to ETS?


title: Separate opening and extension penalties for indels

text: ClustalG, the social science version of ClustalX, can use different 'opening' and 'extension' penalties for indels, so that an indel operation can have more weight when it adds the first element for a new event or gets rid of the last element for an event, instead of having the same weight during the duration of an event. Does TraMineR provide the same way to set up the indel costs?